gestural_alignment_rate
Last updated: 2025-09-05
library(tidyverse)
library(plotrix) #std.error
library(hypr) #for generating contrast coding
library(ggthemes)
library(ggh4x)
library(dagitty) #dagitty for causal diagrams
library(DiagrammeR) #grViz for causal diagrams
library(DiagrammeRsvg)
library(rsvg)
library(plotly) #ggplotly for interactive plots
library(brms) #bayesian regression
library(bayestestR) #p_direction, hdi
library(tidybayes) #add_epred_draws
library(emmeans)
library(doParallel) #set the number of cores
library(ppcor) #partial correlation
# # for power analysis
# library(lme4)
# library(mixedpower)
### Set global options
options(digits = 3) # set the default number of digits to 3
cores = as.integer(detectCores(logical = FALSE) * 0.7) # set the number of cores to use
registerDoParallel(cores=cores) # register the number of cores to use for parallel processing
options(mc.cores = cores)
options(brms.backend = "cmdstanr") #this will speed up the model fitting
### MCMC options
niter = 20000 #number of iterations
nwu = 2000 #number of warmups
### Rmd settings
knitr::opts_chunk$set(echo=TRUE, warning=FALSE, message=FALSE)### Custom functions
########### pp_check_each_round ############
pp_check_each_round <- function(m, data, round_info) {
df_sub = filter(data, round == round_info)
p = pp_check(m,
type = "bars",
ndraws = 100,
newdata = df_sub) +
coord_cartesian(xlim = c(0, 10)) +
ggtitle(round_info)
return(p)
}
########### pp_check_each_condition ############
pp_check_each_condition <- function(m, data, condition_info) {
df_sub = filter(data, condition == condition_info)
p = pp_check(m,
type = "bars",
ndraws = 100,
newdata = df_sub) +
coord_cartesian(xlim = c(0, 10)) +
ggtitle(condition_info)
return(p)
}
########### plot_posterior ############
plot_posterior <- function(model, model2=NA, model3=NA,
interaction=FALSE, include_intercept=FALSE,
xlim_cond=1.5, xlim_round=2, xlim_lex=0.15){
### extract the posterior draws
posterior_beta1 <- model %>%
gather_draws(`b_.*`, regex = TRUE) %>%
mutate(intercept = str_detect(.variable, "Intercept"),
component = ifelse(str_detect(.variable, ":"), "Interaction",
ifelse(str_detect(.variable, "round"), "Round",
ifelse(str_detect(.variable, "Intercept"), "Intercept",
ifelse(str_detect(.variable, "lex_align"), "N. lex align",
"Visibility")))))
if (length(model2) == 1){ #if model2 is NA
posterior_beta = posterior_beta1
} else {
posterior_beta1 <- posterior_beta1 %>%
filter(.variable != "b_num_lex_align")
posterior_beta2 <- model2 %>%
gather_draws(`b_.*`, regex = TRUE) %>%
filter(.variable == "b_num_lex_align") %>%
mutate(component = "N. lex align")
posterior_beta3 <- model3 %>%
gather_draws(`b_.*`, regex = TRUE) %>%
filter(.variable == "b_condition_sumAO_Sym") %>%
mutate(component = "Visibility")
posterior_beta = rbind(posterior_beta1, posterior_beta2, posterior_beta3)
}
posterior_beta = posterior_beta %>%
mutate(.variable = recode(.variable,
"b_Intercept" = "Intercept",
"b_conditionAsymAV" = "SymAV--AsymAV",
"b_conditionAO" = "SymAV--AO",
"b_conditionAsym_Sym" = "AsymAV--SymAV",
"b_conditionAO_Asym" = "AO--AsymAV",
"b_condition_sumAO_Sym" = "AO--SymAV",
"b_num_lex_align" = "N. lex align",
"b_round_c" = "Centered round",
"b_log_round_c" = "Centered log(round)",
"b_roundR12" = "R1--R2",
"b_roundR23" = "R2--R3",
"b_roundR34" = "R3--R4",
"b_roundR45" = "R4--R5",
"b_roundR56" = "R5--R6",
"b_conditionAsymAV:round1" = "SymAV--AsymAV: R1--R2",
"b_conditionAsymAV:round2" = "SymAV--AsymAV: R2--R3",
"b_conditionAsymAV:round3" = "SymAV--AsymAV: R3--R4",
"b_conditionAsymAV:round4" = "SymAV--AsymAV: R4--R5",
"b_conditionAsymAV:round5" = "SymAV--AsymAV: R5--R6",
"b_conditionAO:round1" = "SymAV--AO: R1--R2",
"b_conditionAO:round2" = "SymAV--AO: R2--R3",
"b_conditionAO:round3" = "SymAV--AO: R3--R4",
"b_conditionAO:round4" = "SymAV--AO: R4--R5",
"b_conditionAO:round5" = "SymAV--AO: R5--R6",
"b_conditionAsym_Sym:round_c" = "Centered round: Asym--Sym",
"b_conditionAO_Sym:round_c" = "Centered round: AO--Sym",
"b_conditionAO_Asym:round_c" = "Centered round: AO--Asym",
"b_conditionAsym_Sym:log_round_c" = "Centered log(round): Asym--Sym",
"b_conditionAO_Sym:log_round_c" = "Centered log(round): AO--Sym",
"b_conditionAO_Asym:log_round_c" = "Centered log(round): AO--Asym"),
.variable = factor(.variable,
levels = c("AO--SymAV", "AO--AsymAV", "AsymAV--SymAV",
"R1--R2", "R2--R3", "R3--R4", "R4--R5", "R5--R6",
"N. lex align")),
component = factor(component,
levels = c("Intercept", "Visibility", "Round",
"Interaction", "N. lex align")))
if (include_intercept == F){
posterior_beta = posterior_beta %>% filter(component != "Intercept")
}
### change variables if only main effects are plotted
if (interaction == F) {
posterior_beta = filter(posterior_beta, !str_detect(.variable, ":"))
fill_manual_values = c("steelblue", "steelblue", "steelblue")
} else{
fill_manual_values = c("steelblue", "steelblue", "steelblue", "steelblue")
}
### plot the posterior distributions
p_posterior = ggplot(posterior_beta,
aes(x = .value, y = fct_rev(.variable),
fill = component)) +
geom_vline(xintercept = 0, size = 1) +
stat_halfeye(aes(slab_alpha = intercept),
normalize = "panels",
slab_alpha = .5,
.width = c(0.89, 0.95),
point_interval = "median_hdi") +
scale_fill_manual(values = fill_manual_values) +
scale_slab_alpha_discrete(range = c(0.8, 0.4)) +
guides(fill = "none", slab_alpha = "none") +
labs(x = "Coefficient", y = "Effect") +
theme_clean(base_size = 15) +
theme(axis.text.x = element_text(colour = "black", size = 14),
axis.text.y = element_text(colour = "black", size = 14),
axis.title = element_text(size = 15, face = 'bold'),
axis.title.x = element_text(vjust = -2),
axis.title.y = element_text(vjust = 2),
legend.position = "none",
strip.text = element_text(size = 15, face = 'bold'),
strip.background = element_blank(),
panel.grid.major.x = element_line(color = "grey90",
linetype = "solid",
size = 0.5),
panel.grid.major.y = element_line(color = "grey90",
linetype = "solid",
size = 0.5),
plot.background = element_blank(),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines")) +
facet_wrap(vars(component), ncol = 3, scales = "free") +
facetted_pos_scales(
x = list(
scale_x_continuous(limits = c(-xlim_cond, xlim_cond)),
scale_x_continuous(limits = c(-xlim_round, xlim_round)),
scale_x_continuous(limits = c(-xlim_lex, xlim_lex),
breaks = c(-0.1, 0, 0.1))))
return(p_posterior)
}
########### pp_update_plot ############
pp_update_plot <- function(post_sample, model_type="nb", interaction=TRUE){
sum = ifelse("b_condition_sumAO_Sym" %in% colnames(post_sample), T, F)
round_bd = ifelse("b_roundR12" %in% colnames(post_sample), T, F)
intercept = ggplot(post_sample) +
geom_density(aes(prior_Intercept), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_Intercept), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Intercept') +
theme_classic()
### Visibility condition
if (sum == F){
cond1 = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_conditionAsym_Sym), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Asym--Sym') +
theme_classic()
cond2 = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_conditionAO_Asym), fill="#FC4E07", color="black",alpha=0.6) +
xlab('AO--Asym') +
theme_classic()
} else {
cond1 = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_condition_sumAO_Sym), fill="#FC4E07", color="black",alpha=0.6) +
xlab('AO--Sym') +
theme_classic()
cond2 = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_condition_sumAsym_Sym), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Asym--Sym') +
theme_classic()
}
### Round
if (interaction) {
if (round_bd){
r1 = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_roundR12), fill="#FC4E07", color="black",alpha=0.6) +
xlab('R1--R2') +
theme_classic()
r2 = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_roundR23), fill="#FC4E07", color="black",alpha=0.6) +
xlab('R2--R3') +
theme_classic()
r3 = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_roundR34), fill="#FC4E07", color="black",alpha=0.6) +
xlab('R3--R4') +
theme_classic()
r4 = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_roundR45), fill="#FC4E07", color="black",alpha=0.6) +
xlab('R4--R5') +
theme_classic()
r5 = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_roundR56), fill="#FC4E07", color="black",alpha=0.6) +
xlab('R5--R6') +
theme_classic()
} else {
if ("b_round_c" %in% colnames(post_sample)) {
round = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_round_c), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Centered round') +
theme_classic()
if (sum == F){
cond1_round = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(`b_conditionAsym_Sym:round_c`), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Centered Round: Asym--Sym') +
theme_classic()
cond2_round = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(`b_conditionAO_Asym:round_c`), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Centered Round: AO--Asym') +
theme_classic()
} else {
cond1_round = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(`b_condition_sumAO_Sym:round_c`), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Centered Round: AO--Sym') +
theme_classic()
cond2_round = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(`b_condition_sumAsym_Sym:round_c`), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Centered Round: Asym--Sym') +
theme_classic()
}
} else {
round = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_log_round_c), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Centered log(round)') +
theme_classic()
if (sum == F){
cond1_round = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(`b_conditionAsym_Sym:log_round_c`), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Centered log(round): Asym--Sym') +
theme_classic()
cond2_round = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(`b_conditionAO_Asym:log_round_c`), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Centered log(round): AO--Asym') +
theme_classic()
} else {
cond1_round = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(`b_condition_sumAO_Sym:log_round_c`), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Centered log(round): AO--Sym') +
theme_classic()
cond2_round = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(`b_condition_sumAsym_Sym:log_round_c`), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Centered log(round): Asym--Sym') +
theme_classic()
}
}
}
}
if (model_type == "nb"){
shape = ggplot(post_sample) +
geom_density(aes(prior_shape), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(shape), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Shape') +
scale_x_continuous(limits = c(0, 10)) +
theme_classic()}
else if (model_type == "zinb") {
shape = ggplot(post_sample) +
geom_density(aes(prior_shape), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(shape), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Shape') +
scale_x_continuous(limits = c(0, 10)) +
theme_classic()
zi = ggplot(post_sample) +
geom_density(aes(prior_zi), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(zi), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Zero-inflation') +
theme_classic()}
else if (model_type == "zibt"){
phi = ggplot(post_sample) +
geom_density(aes(prior_phi), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(phi), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Precision') +
theme_classic()
zoi = ggplot(post_sample) +
geom_density(aes(prior_zoi), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(zoi), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Zero-inflation') +
theme_classic()
coi = ggplot(post_sample) +
geom_density(aes(prior_coi), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(coi), fill="#FC4E07", color="black",alpha=0.6) +
xlab('One-inflation') +
theme_classic()
}
### display the plots
if (interaction==F){
if (model_type=="nb"){
gridExtra::grid.arrange(intercept, cond1, cond2, shape, ncol=2)}
else if (model_type=="zinb"){
gridExtra::grid.arrange(intercept, cond1, cond2, shape, zi, ncol=3)}
else if (model_type=="zibt"){
gridExtra::grid.arrange(intercept, cond1, cond2, phi, zoi, coi, ncol=3)
}
} else {
if (round_bd == T){
if (model_type=="nb"){
gridExtra::grid.arrange(intercept, cond1, cond2, shape, ncol=3)}
else if (model_type=="zinb"){
gridExtra::grid.arrange(intercept, cond1, cond2, shape, zi, ncol=3)}
else if (model_type=="zibt"){
gridExtra::grid.arrange(intercept, cond1, cond2, phi, zoi, coi, ncol=3)
}
} else {
if (model_type=="nb"){
gridExtra::grid.arrange(intercept, cond1, cond2, round,
cond1_round, cond2_round, shape, ncol=3)}
else if (model_type=="zinb"){
gridExtra::grid.arrange(intercept, cond1, cond2, round,
cond1_round, cond2_round, shape, zi, ncol=3)}
else if (model_type=="zibt"){
gridExtra::grid.arrange(intercept, cond1, cond2, round,
cond1_round, cond2_round, phi, zoi, coi, ncol=3)
}
}
}
}1 =====Data preparation=====
1.1 Load data
### trial_info.csv
df_trial_info = read.csv("data/trial_info.csv", stringsAsFactors = TRUE) %>%
filter(num_words != 0) %>% # remove trials that were accidentally skipped) %>%
mutate(pair = factor(pair),
round_c = as.integer(round) - mean(as.integer(round)),
round_n = factor(round),
round = factor(paste0("R", round)),
log_round = log(as.integer(round)),
log_round_c = log_round - mean(log_round),
condition = factor(condition,
levels = c("Sym", "Asym", "AO"),
labels = c("SymAV", "AsymAV", "AO")),
condition_sum = condition,
duration_s = duration_ms / 1000,
n_words_100 = num_words / 100,
n_iconic_per_100words = num_iconic_gestures / n_words_100,
n_iconic_100 = num_iconic_gestures / 100) %>%
dplyr::select(pair, condition, round, round_n,
round_c, log_round, log_round_c,
duration_ms, duration_s, duration_m,
director:n_iconic_100) %>%
rename(trial_duration_s = duration_s, trial_duration_m = duration_m)
head(df_trial_info)### condition info
df_condition_info = read.csv("data/condition_info.csv", stringsAsFactors = TRUE) %>%
mutate(pair = factor(pair),
condition = factor(condition,
levels = c("Sym", "Asym", "AO"),
labels = c("SymAV", "AsymAV", "AO")))2 Demographics
2.1 Participant
### participants
df_participant = read_tsv("data/participant_info.txt") %>%
mutate(pair = factor(as.integer(pair)),
condition = factor(condition))
df_participant_used = df_participant %>%
filter(used == "Y") %>%
mutate(condition = factor(condition, levels = c("Sym", "Asym", "AO")))
df_participant_usedmean_age = mean(df_participant_used$age)
sd_age = sd(df_participant_used$age)
range_age = paste(range(df_participant_used$age)[1], range(df_participant_used$age)[2], sep = " - ")
n_participants = nrow(df_participant_used)
n_female = nrow(filter(df_participant_used, gender == "F"))
n_male = nrow(filter(df_participant_used, gender == "M"))
data.frame(mean_age, sd_age, range_age, n_participants, n_female, n_male)### dyads
# make a dataframe with dyad information (for each pair, gender & age of speaker 1 & 2)
df_dyad = df_participant_used %>%
dplyr::select(-notes) %>%
pivot_wider(names_from = speaker,
values_from = c("gender", "age")) %>%
mutate(mixed_gender = ifelse(gender_A == gender_B, "N", "Y"),
mixed_gender = ifelse(is.na(mixed_gender), "Y", mixed_gender),
gender_dyad = paste0(gender_A, "_", gender_B))
df_dyad %>%
group_by(condition, mixed_gender) %>%
summarize(n = n())df_dyad %>%
group_by(condition, mixed_gender, gender_dyad) %>%
summarize(n = n())2.2 Trial information & iconic gestures
summarize_demographics <- function(df) {
df %>%
summarize(total_s = sum(trial_duration_s),
total_m = total_s / 60,
### trial duration ###
mean_trial_dur_s = mean(trial_duration_s),
mean_trial_dur_m = mean(trial_duration_m),
sd_trial_dur_m = sd(trial_duration_m),
lci_trial_dur_m = mean_trial_dur_m - 1.96 * sd_trial_dur_m / sqrt(n()),
uci_trial_dur_m = mean_trial_dur_m + 1.96 * sd_trial_dur_m / sqrt(n()),
### words ###
# number of words
num_words_total = sum(num_words),
mean_num_words = mean(num_words),
num_words_100 = mean(num_words / 100),
sd_num_words = sd(num_words),
lci_num_words = mean_num_words - 1.96 * sd_num_words / sqrt(n()),
uci_num_words = mean_num_words + 1.96 * sd_num_words / sqrt(n()),
num_words_per_min = num_words_total / total_m,
# number of content words
num_content_total = sum(num_content_words),
mean_num_content = mean(num_content_words),
sd_num_content = sd(num_content_words),
lci_num_content = mean_num_content - 1.96 * sd_num_content / sqrt(n()),
uci_num_content = mean_num_content + 1.96 * sd_num_content / sqrt(n()),
num_content_per_min = num_content_total / total_m,
### lexical alignment ###
# raw frequency
num_lex_align_total = sum(num_lex_align, na.rm = T),
mean_num_lex_align = mean(num_lex_align, na.rm = T),
sd_num_lex_align = sd(num_lex_align, na.rm = T),
lci_num_lex_align = mean_num_lex_align - 1.96 * sd_num_lex_align / sqrt(n()),
uci_num_lex_align = mean_num_lex_align + 1.96 * sd_num_lex_align / sqrt(n()),
num_lex_align_per_min = num_lex_align_total / total_m,
# rate per 100 words
mean_lex_align_per_100words = mean(num_lex_align / num_words_100, na.rm=T),
sd_lex_align_per_100words = sd(num_lex_align / num_words_100, na.rm=T),
se_lex_align_per_100words = std.error(num_lex_align / num_words_100, na.rm=T),
lci_lex_align_per_100words = mean_lex_align_per_100words - 1.96 * se_lex_align_per_100words,
uci_lex_align_per_100words = mean_lex_align_per_100words + 1.96 * se_lex_align_per_100words,
### iconic gestures ###
# raw frequency
num_iconic_total = sum(num_iconic_gestures, na.rm = T),
mean_num_iconic = mean(num_iconic_gestures, na.rm = T),
sd_num_iconic = sd(num_iconic_gestures, na.rm = T),
lci_num_iconic = mean_num_iconic - 1.96 * sd_num_iconic / sqrt(n()),
uci_num_iconic = mean_num_iconic + 1.96 * sd_num_iconic / sqrt(n()),
num_iconic_per_min = num_iconic_total / total_m,
# rate per 100 words
mean_iconic_per_100words = mean(num_iconic_gestures / num_words_100, na.rm=T),
sd_iconic_per_100words = sd(num_iconic_gestures / num_words_100, na.rm=T),
se_iconic_per_100words = std.error(num_iconic_gestures / num_words_100, na.rm=T),
lci_iconic_per_100words = mean_iconic_per_100words - 1.96 * se_iconic_per_100words,
uci_iconic_per_100words = mean_iconic_per_100words + 1.96 * se_iconic_per_100words,
### gestural alignment ###
# raw frequency
num_gest_align_total = sum(num_gestural_align, na.rm = T),
mean_num_gest_align = mean(num_gestural_align, na.rm = T),
sd_num_gest_align = sd(num_gestural_align, na.rm = T),
lci_num_gest_align = mean_num_gest_align - 1.96 * sd_num_gest_align / sqrt(n()),
uci_num_gest_align = mean_num_gest_align + 1.96 * sd_num_gest_align / sqrt(n()),
num_gest_align_per_min = num_gest_align_total / total_m,
# rate per 100 words
mean_gest_align_per_100words = mean(num_gestural_align / num_words_100, na.rm=T),
sd_gest_align_per_100words = sd(num_gestural_align / num_words_100, na.rm=T),
se_gest_align_per_100words = std.error(num_gestural_align / num_words_100, na.rm=T),
lci_gest_align_per_100words = mean_gest_align_per_100words - 1.96 * se_gest_align_per_100words,
uci_gest_align_per_100words = mean_gest_align_per_100words + 1.96 * se_gest_align_per_100words,
# rate per iconic gestures
mean_gest_align_prop = mean(num_gestural_align / num_iconic_gestures, na.rm=T),
sd_gest_align_per_prop = sd(num_gestural_align / num_iconic_gestures, na.rm=T),
se_gest_align_per_prop = std.error(num_gestural_align / num_iconic_gestures, na.rm=T),
lci_gest_align_per_prop = mean_gest_align_prop - 1.96 * se_gest_align_per_prop,
uci_gest_align_per_prop = mean_gest_align_prop + 1.96 * se_gest_align_per_prop,
### number of trials ###
trial_n = n()) %>%
ungroup()
}
summarize_dur <- function(df){
df %>%
summarize(mean_dur_m = mean(total_m),
sd_dur_m = sd(total_m),
se_dur_m = std.error(total_m),
lci_dur_m = mean_dur_m - 1.96 * se_dur_m,
uci_dur_m = mean_dur_m + 1.96 * se_dur_m) %>%
ungroup()
}2.2.1 mean by condition
### demographics by pair
trial_info_pair = df_trial_info %>%
group_by(pair) %>%
summarize_demographics() %>%
left_join(., df_condition_info) %>%
dplyr::select(pair, condition, total_s:trial_n)
### calculate mean experiment duration
mean_dur_cond = trial_info_pair %>%
group_by(condition) %>%
summarize_dur()
### summary statistics
trial_info_cond = df_trial_info %>%
group_by(condition) %>%
summarize_demographics() %>%
left_join(., mean_dur_cond) %>%
dplyr::select(condition, total_m, mean_dur_m:uci_dur_m,
everything(),
-ends_with("_s"))
trial_info_cond2.2.2 mean by round
trial_info_pair_round = df_trial_info %>%
group_by(pair, round, round_n) %>%
summarize_demographics() %>%
left_join(., df_condition_info)
### calculate mean round duration
mean_dur_round = trial_info_pair_round %>%
group_by(round, round_n) %>%
summarize_dur()
trial_info_round = df_trial_info %>%
group_by(round, round_n) %>%
summarize_demographics() %>%
left_join(., mean_dur_round) %>%
dplyr::select(round, round_n, total_m, mean_dur_m:uci_dur_m,
everything(),
-ends_with("_s"))
trial_info_round2.2.3 mean by condition x round
### calculate mean duration by condition x round
mean_dur_cond_round = trial_info_pair_round %>%
group_by(condition, round, round_n) %>%
summarize_dur()
trial_info_cond_round = df_trial_info %>%
group_by(condition, round, round_n) %>%
summarize_demographics() %>%
left_join(., mean_dur_cond_round) %>%
dplyr::select(condition, round, round_n,
total_m, mean_dur_m:uci_dur_m,
everything(),
-ends_with("_s"))
trial_info_cond_round3 =====Contrast coding=====
### visibility condition: difference coding
h_cond = hypr(AO_Asym = AsymAV ~ AO,
Asym_Sym = SymAV ~ AsymAV,
levels = levels(df_trial_info$condition))
h_cond## hypr object containing 2 null hypotheses:
## H0.AO_Asym: 0 = AsymAV - AO
## H0.Asym_Sym: 0 = SymAV - AsymAV
##
## Call:
## hypr(AO_Asym = ~AsymAV - AO, Asym_Sym = ~SymAV - AsymAV, levels = c("SymAV",
## "AsymAV", "AO"))
##
## Hypothesis matrix (transposed):
## AO_Asym Asym_Sym
## SymAV 0 1
## AsymAV 1 -1
## AO -1 0
##
## Contrast matrix:
## AO_Asym Asym_Sym
## SymAV 1/3 2/3
## AsymAV 1/3 -1/3
## AO -2/3 -1/3contrasts(df_trial_info$condition) = contr.hypothesis(h_cond)
### visibility condition: sum coding
h_cond = hypr(AO_Sym = SymAV ~ AO,
Asym_Sym = SymAV ~ AsymAV,
levels = levels(df_trial_info$condition))
h_cond## hypr object containing 2 null hypotheses:
## H0.AO_Sym: 0 = SymAV - AO
## H0.Asym_Sym: 0 = SymAV - AsymAV
##
## Call:
## hypr(AO_Sym = ~SymAV - AO, Asym_Sym = ~SymAV - AsymAV, levels = c("SymAV",
## "AsymAV", "AO"))
##
## Hypothesis matrix (transposed):
## AO_Sym Asym_Sym
## SymAV 1 1
## AsymAV 0 -1
## AO -1 0
##
## Contrast matrix:
## AO_Sym Asym_Sym
## SymAV 1/3 1/3
## AsymAV 1/3 -2/3
## AO -2/3 1/3contrasts(df_trial_info$condition_sum) = contr.hypothesis(h_cond)
### round
bacward_diff = hypr(R12 = R2 ~ R1,
R23 = R3 ~ R2,
R34 = R4 ~ R3,
R45 = R5 ~ R4,
R56 = R6 ~ R5,
levels = levels(df_trial_info$round))
bacward_diff## hypr object containing 5 null hypotheses:
## H0.R12: 0 = R2 - R1
## H0.R23: 0 = R3 - R2
## H0.R34: 0 = R4 - R3
## H0.R45: 0 = R5 - R4
## H0.R56: 0 = R6 - R5
##
## Call:
## hypr(R12 = ~R2 - R1, R23 = ~R3 - R2, R34 = ~R4 - R3, R45 = ~R5 -
## R4, R56 = ~R6 - R5, levels = c("R1", "R2", "R3", "R4", "R5",
## "R6"))
##
## Hypothesis matrix (transposed):
## R12 R23 R34 R45 R56
## R1 -1 0 0 0 0
## R2 1 -1 0 0 0
## R3 0 1 -1 0 0
## R4 0 0 1 -1 0
## R5 0 0 0 1 -1
## R6 0 0 0 0 1
##
## Contrast matrix:
## R12 R23 R34 R45 R56
## R1 -5/6 -2/3 -1/2 -1/3 -1/6
## R2 1/6 -2/3 -1/2 -1/3 -1/6
## R3 1/6 1/3 -1/2 -1/3 -1/6
## R4 1/6 1/3 1/2 -1/3 -1/6
## R5 1/6 1/3 1/2 2/3 -1/6
## R6 1/6 1/3 1/2 2/3 5/6contrasts(df_trial_info$round) = contr.hypothesis(bacward_diff)4 =====Number of iconic gestures=====
4.1 —DataViz: frequency—
4.1.1 bp: total by condition
bp_iconic_by_cond = ggplot(data=trial_info_pair,
aes(x=condition, y=num_iconic_total, fill=condition)) +
geom_jitter(aes(color = pair),
size = 1, alpha = 1,
width = 0.07, height = 0) +
geom_boxplot(width = .2,
outlier.shape = NA, alpha = 0.7) +
scale_fill_manual(values = c("#ED6B06", "#00786A", "darkgrey")) +
labs(x = "Visibility",
y = "Total N of iconic gestures per pair") +
theme_clean(base_size = 15) +
theme(axis.text.x = element_text(colour = "black", size = 14),
axis.text.y = element_text(colour = "black", size = 14),
axis.title = element_text(size = 15, face = 'bold'),
axis.title.x = element_text(vjust = -2),
axis.title.y = element_text(vjust = 2),
legend.position = "none",
strip.text = element_text(size = 15, face = 'bold'),
plot.background = element_blank(),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines"))
ggplotly(bp_iconic_by_cond)4.1.2 bp: mean by condition
bp_mean_iconic_by_cond = ggplot(data=trial_info_pair,
aes(x=condition, y=mean_num_iconic, fill=condition)) +
geom_jitter(aes(color = pair),
size = 1, alpha = 1,
width = 0.07, height = 0) +
geom_boxplot(width = .2,
outlier.shape = NA, alpha = 0.7) +
geom_point(data = trial_info_cond,
aes(y = mean_num_iconic),
size = 2, shape = 4) +
scale_y_continuous(limits = c(0, 4), breaks = seq(0, 4, 1)) +
scale_fill_manual(values = c("#ED6B06", "#00786A", "darkgrey")) +
labs(x = "Visibility",
y = "Mean N of iconic gestures per trial") +
theme_clean(base_size = 15) +
theme(axis.text.x = element_text(colour = "black", size = 14),
axis.text.y = element_text(colour = "black", size = 14),
axis.title = element_text(size = 15, face = 'bold'),
axis.title.x = element_text(vjust = -2),
axis.title.y = element_text(vjust = 2),
legend.position = "none",
strip.text = element_text(size = 15, face = 'bold'),
plot.background = element_blank(),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines"))
ggplotly(bp_mean_iconic_by_cond)4.1.3 bp: mean by condition x round
bp_mean_iconic_by_round_cond = ggplot(data=df_trial_info,
aes(x=round, y=num_iconic_gestures, fill=condition)) +
geom_boxplot(outlier.shape = NA,
alpha = 0.7) +
scale_fill_manual(values = c("#ED6B06", "#00786A", "darkgrey")) +
scale_y_continuous(limits = c(0, 15), breaks = seq(0, 15, 5)) +
labs(x = "Round",
y = "Mean N of iconic gestures per trial") +
theme_clean(base_size = 15) +
theme(axis.text.x = element_text(colour = "black", size = 14),
axis.text.y = element_text(colour = "black", size = 14),
axis.title = element_text(size = 15, face = 'bold'),
axis.title.x = element_text(vjust = -2),
axis.title.y = element_text(vjust = 2),
legend.position = "top",
strip.text = element_text(size = 15, face = 'bold'),
plot.background = element_blank(),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines"))
bp_mean_iconic_by_round_cond
The figure shows that the decline in iconic gestures across rounds does not follow a linear pattern. This suggests that log-transformed round may improve the model fit.
4.2 —Zero-inflated Poisson regression models—
4.2.1 Prior specification
We will use weakly informative priors for the regression each parameter. To specify priors for the intercept, which represents the expected number of iconic gestures in the SymAV condition, we will refer to the number of iconic gestures reported in Akamine et al. (2024). In the paper, the authors analyzed data from 19 dyads involving in the same task as the current study but in co-present interaction. The total number of iconic gestures reported was 4413, which was collected from 19 dyads, each performing 96 trials. Therefore, the expected number of iconic gestures per trial is \(4413 / (19 * 96) = 2.42\).
As Poisson regression models use a log-link, we need to specify priors on the log scale. Taking the log of 2.42, we get 0.88. In order to allow data to inform the posterior, we will use a normal prior with a mean of 0.88 and a relatively wide standard deviation of 3.
For the fixed effects, we will set unbiased priors with a mean of 0 and a standard deviation of 2.
For the standard deviation of the random effects, we set the prior to be normal with mean 0 and standard deviation 2. For the correlation between the random effects, we set the prior to be LKJ(2).
For the models including the round as fixed effects (whether backward-difference coded or centered + log-transformed), the intercept will represent the mean expected number of iconic gestures (ground mean) in the SymAV condition. As the meaning of the intercept doesn’t change when adding the round variable, we use the same prior for the intercept.
priors_rint = c(
prior(normal(0.88, 2), class = Intercept),
prior(normal(0, 2), class = b),
prior(normal(0, 2), class = sd))
priors_rslope = c(
prior(normal(0.88, 2), class = Intercept),
prior(normal(0, 2), class = b),
prior(normal(0, 2), class = sd),
prior(lkj(2), class = cor))4.2.2 Model 1: condition (pair)
First, we will model the number of iconic gestures per trial with fixed effects for condition. As we expect no gestures in many trials based on previous study (Akamine et al., 2024), we will use a zero-inflated Poisson regression model with a log link function.
4.2.2.1 Prior predictive check
mp_iconic_cond_prior = brm(num_iconic_gestures ~ 1 + condition +
(1|pair) + (1|target),
family = poisson(),
prior = priors_rint,
sample_prior = "only",
data = df_trial_info,
file = "models/mp_iconic_cond_prior")
pp_check(mp_iconic_cond_prior, ndraws = 100, type = "bars") +
coord_cartesian(xlim = c(0, 20),
ylim = c(0, 4000))
The prior predictive check shows that the model generates data that are similar to the observed data. The model generates a far wider range of data than the observed data, allowing the model to be informed by the data.
Next, we will fit the model to the observed data.
4.2.2.2 Fit the model
mp_iconic_cond = brm(num_iconic_gestures ~ 1 + condition +
(1|pair) + (1|target),
family = zero_inflated_poisson(),
prior = priors_rint,
data = df_trial_info,
sample_prior = T,
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/mp_iconic_cond")
summary(mp_iconic_cond)## Family: zero_inflated_poisson
## Links: mu = log; zi = identity
## Formula: num_iconic_gestures ~ 1 + condition + (1 | pair) + (1 | target)
## Data: df_trial_info (Number of observations: 4316)
## Draws: 4 chains, each with iter = 20000; warmup = 2000; thin = 1;
## total post-warmup draws = 72000
##
## Multilevel Hyperparameters:
## ~pair (Number of levels: 45)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.55 0.07 0.44 0.69 1.00 10772 22264
##
## ~target (Number of levels: 16)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.21 0.05 0.14 0.33 1.00 12904 27375
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.91 0.10 0.71 1.11 1.00 8817 16940
## conditionAO_Asym 0.29 0.21 -0.12 0.70 1.00 8349 15439
## conditionAsym_Sym -0.07 0.20 -0.48 0.33 1.00 8526 15481
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## zi 0.41 0.01 0.39 0.43 1.00 56110 51525
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).bayestestR::hdi(mp_iconic_cond)bayestestR::hdi(mp_iconic_cond, ci = 0.89)4.2.2.3 Posterior predictive check
pp_check_sym = pp_check_each_condition(mp_iconic_cond, df_trial_info, "SymAV")
pp_check_asym = pp_check_each_condition(mp_iconic_cond, df_trial_info, "AsymAV")
pp_check_ao = pp_check_each_condition(mp_iconic_cond, df_trial_info, "AO")
gridExtra::grid.arrange(pp_check_sym, pp_check_asym, pp_check_ao, ncol = 2)
The posterior predictive check shows that the model generates data that are similar to the observed data. This suggests that the model is a good fit to the data.
4.2.2.4 Model Effect
Here, we will compare which regression model is suitable for the data, poisson or negative binomial. We will follow the procedure introduced in Winter & Burkner (2021).
First thing first, we will fit the negative binomial regression model to the data. Then we compare the predictive power of the models using loo_compare() function.
mnb_iconic_cond = brm(num_iconic_gestures ~ 1 + condition +
(1|pair) + (1|target),
family = negbinomial(),
prior = priors_rint,
data = df_trial_info,
sample_prior = T,
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/mnb_iconic_cond")
mnb_iconic_cond_zi = brm(num_iconic_gestures ~ 1 + condition +
(1|pair) + (1|target),
family = zero_inflated_negbinomial(),
prior = priors_rint,
data = df_trial_info,
sample_prior = T,
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/mnb_iconic_cond_zi")
### loo compare
if (!file.exists("models/loo_pois_iconic.rds")){
pois_loo = loo(mp_iconic_cond)
saveRDS(pois_loo, file = "models/loo_pois_iconic.rds")
}
if (!file.exists("models/loo_nb_iconic.rds")){
nb_loo = loo(mnb_iconic_cond)
saveRDS(nb_loo, file = "models/loo_nb_iconic.rds")
}
if (!file.exists("models/loo_zinb_iconic.rds")){
zinb_loo = loo(mnb_iconic_cond_zi)
saveRDS(zinb_loo, file = "models/loo_zinb_iconic.rds")
}
pois_loo = readRDS("models/loo_pois_iconic.rds")
nb_loo = readRDS("models/loo_nb_iconic.rds")
zinb_loo = readRDS("models/loo_zinb_iconic.rds")
loo_compare(pois_loo, nb_loo, zinb_loo)## elpd_diff se_diff
## mnb_iconic_cond 0.0 0.0
## mnb_iconic_cond_zi -0.9 1.5
## mp_iconic_cond -947.7 102.0The loo_compare() function shows that the negative binomial regression has a larger predictive power (elpd_diff) and smaller standard error (se_diff) than the zero-inflated negative binomial regression and the Poisson regression. Therefore, we will use the negative binomial regression model for further analyses.
4.3 —Negative binomial regression models—
4.3.1 Model 2: condition x round (pair)
Next, we will run the negative binomial regression model to predict the number of iconic gestures per trial by condition and round. We will use the priors_round for the model.
4.3.1.1 Model effect (log round)
We saw earlier that the decline in # iconic gestures doesn’t follow a linear pattern: The decline is stepper in the initial rounds and then plateaus afterwards. To capture this trend better without making 5 levels for Round (which is the case for bd coding), we will check if the model performance is comparable for backward-difference coded round and centered log-transformed round.
mnb_iconic_cond_round = brm(num_iconic_gestures ~ 1 + condition * round +
(1+round|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope,
data = df_trial_info,
sample_prior = T,
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/mnb_iconic_cond_round")
mnb_iconic_cond_round_c = brm(num_iconic_gestures ~ 1 + condition * round_c +
(1+round_c|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope,
data = df_trial_info,
sample_prior = T,
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/mnb_iconic_cond_round_c")
mnb_iconic_cond_round_log = brm(num_iconic_gestures ~ 1 + condition * log_round_c +
(1+log_round_c|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope,
data = df_trial_info,
sample_prior = T,
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/mnb_iconic_cond_round_log")
### loo compare
if (!file.exists("models/loo_nb_iconic_cond_round.rds")){
nb_cond_round_loo = loo(mnb_iconic_cond_round)
saveRDS(nb_cond_round_loo, file = "models/loo_nb_iconic_cond_round.rds")
}
if (!file.exists("models/loo_nb_iconic_cond_round_c.rds")){
nb_cond_round_c_loo = loo(mnb_iconic_cond_round_c)
saveRDS(nb_cond_round_c_loo, file = "models/loo_nb_iconic_cond_round_c.rds")
}
if (!file.exists("models/loo_nb_iconic_cond_round_log.rds")){
nb_cond_round_log_loo = loo(mnb_iconic_cond_round_log)
saveRDS(nb_cond_round_log_loo, file = "models/loo_nb_iconic_cond_round_log.rds")
}
nb_cond_round_loo = readRDS("models/loo_nb_iconic_cond_round.rds")
nb_cond_round_c_loo = readRDS("models/loo_nb_iconic_cond_round_c.rds")
nb_cond_round_log_loo = readRDS("models/loo_nb_iconic_cond_round_log.rds")
loo_compare(nb_cond_round_loo, nb_cond_round_c_loo, nb_cond_round_log_loo)## elpd_diff se_diff
## mnb_iconic_cond_round_log 0.0 0.0
## mnb_iconic_cond_round -13.9 5.3
## mnb_iconic_cond_round_c -38.2 12.8The leave-one-out (LOO) Effect shows that centered log-transformed round provides a larger predictive power (elpd_diff) and smaller standard error (se_diff) compared to the backward-difference coded round or centered round. This is in line with our observation that the decrease in the number of gestures does not follow a linear pattern, which explains why the centered round is disfavored here. Therefore, we will use the centered log round.
4.3.1.2 Prior predictive check
mnb_iconic_cond_round_log = brm(num_iconic_gestures ~ 1 + condition * log_round_c +
(1+log_round_c|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope,
sample_prior = "only",
data = df_trial_info,
file = "models/mnb_iconic_cond_round_log_prior")
pp_check(mnb_iconic_cond_round_log, ndraws = 100, type = "bars") +
coord_cartesian(xlim = c(0, 20),
ylim = c(0, 4000))
The prior predictive check shows that the model generates data that are similar to the observed data.
4.3.1.3 Fit the model
mnb_iconic_cond_round_log = brm(num_iconic_gestures ~ 1 + condition * log_round_c +
(1+log_round_c|pair) + (1|target),
family = negbinomial(),
prior = priors_round,
data = df_trial_info,
sample_prior = T,
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/mnb_iconic_cond_round_log")
summary(mnb_iconic_cond_round_log)## Family: negbinomial
## Links: mu = log; shape = identity
## Formula: num_iconic_gestures ~ 1 + condition * log_round_c + (1 + log_round_c | pair) + (1 | target)
## Data: df_trial_info (Number of observations: 4316)
## Draws: 4 chains, each with iter = 20000; warmup = 2000; thin = 1;
## total post-warmup draws = 72000
##
## Multilevel Hyperparameters:
## ~pair (Number of levels: 45)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.85 0.10 0.68 1.07 1.00 15511
## sd(log_round_c) 0.41 0.06 0.31 0.54 1.00 24199
## cor(Intercept,log_round_c) 0.64 0.11 0.38 0.82 1.00 30037
## Tail_ESS
## sd(Intercept) 26870
## sd(log_round_c) 37632
## cor(Intercept,log_round_c) 44753
##
## ~target (Number of levels: 16)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.30 0.07 0.20 0.45 1.00 18463 33962
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept -0.05 0.15 -0.34 0.25 1.00
## conditionAO_Asym 0.37 0.31 -0.24 0.98 1.00
## conditionAsym_Sym -0.16 0.31 -0.77 0.45 1.00
## log_round_c -1.37 0.07 -1.51 -1.23 1.00
## conditionAO_Asym:log_round_c 0.00 0.17 -0.33 0.34 1.00
## conditionAsym_Sym:log_round_c -0.10 0.17 -0.43 0.23 1.00
## Bulk_ESS Tail_ESS
## Intercept 12579 22655
## conditionAO_Asym 11906 21005
## conditionAsym_Sym 11866 21349
## log_round_c 18656 32465
## conditionAO_Asym:log_round_c 19212 31521
## conditionAsym_Sym:log_round_c 19007 32931
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 2.06 0.13 1.83 2.32 1.00 81667 50845
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).bayestestR::hdi(mnb_iconic_cond_round_log)# bayestestR::hdi(mnb_iconic_cond_round_log, ci = 0.89)Based on the model estimates, the condition does not have a significant effect on the number of iconic gestures produced per trial. However, the round has a significant effect on the number of iconic gestures produced per trial: the number of iconic gestures significantly decreases as the round progresses.
4.3.1.4 Visualize the posterior distributions
# p = plot_posterior(mnb_iconic_cond_round_log, interaction=T, xlim=1)
# p4.3.1.5 Probability of direction
p_direction(mnb_iconic_cond_round_log)The probability of direction demonstrated significant effects for the round only.
4.3.1.6 Prior-posterior update plot
### condition
#Sample the parameters of interest:
posterior <- as_draws_df(mnb_iconic_cond_round_log)
#Plot the prior-posterior update plot for the estimated parameters:
pp_update_intercept = ggplot(posterior) +
geom_density(aes(prior_Intercept), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_Intercept), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Intercept') +
theme_classic()
pp_update_b1 = ggplot(posterior) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_conditionAO_Asym), fill="#FC4E07", color="black",alpha=0.6) +
xlab('AO --Asym') +
theme_classic()
pp_update_b2 = ggplot(posterior) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_conditionAsym_Sym), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Asym--Sym') +
theme_classic()
### round
pp_update_round = ggplot(posterior) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_log_round_c), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Centered log(round)') +
theme_classic()
gridExtra::grid.arrange(pp_update_intercept, pp_update_b1, pp_update_b2, pp_update_round)
4.3.1.7 Posterior predictive check
pp_check_sym = pp_check_each_condition(mnb_iconic_cond_round_log, df_trial_info, "SymAV")
pp_check_asym = pp_check_each_condition(mnb_iconic_cond_round_log, df_trial_info, "AsymAV")
pp_check_ao = pp_check_each_condition(mnb_iconic_cond_round_log, df_trial_info, "AO")
gridExtra::grid.arrange(pp_check_sym, pp_check_asym, pp_check_ao, ncol = 2)
pp_check_r1 = pp_check_each_round(mnb_iconic_cond_round_log, df_trial_info, "R1")
pp_check_r2 = pp_check_each_round(mnb_iconic_cond_round_log, df_trial_info, "R2")
pp_check_r3 = pp_check_each_round(mnb_iconic_cond_round_log, df_trial_info, "R3")
pp_check_r4 = pp_check_each_round(mnb_iconic_cond_round_log, df_trial_info, "R4")
pp_check_r5 = pp_check_each_round(mnb_iconic_cond_round_log, df_trial_info, "R5")
pp_check_r6 = pp_check_each_round(mnb_iconic_cond_round_log, df_trial_info, "R6")
gridExtra::grid.arrange(pp_check_r1, pp_check_r2, pp_check_r3, pp_check_r4, pp_check_r5, pp_check_r6, ncol = 3)
4.3.1.8 Pair-wise Effect
emmeans(mnb_iconic_cond_round_log, pairwise ~ condition)$contrasts## contrast estimate lower.HPD upper.HPD
## SymAV - AsymAV -0.160 -0.769 0.454
## SymAV - AO 0.214 -0.402 0.830
## AsymAV - AO 0.375 -0.240 0.981
##
## Point estimate displayed: median
## Results are given on the log (not the response) scale.
## HPD interval probability: 0.95# emmeans(mnb_iconic_cond_round_log, pairwise ~ condition, level = 0.89)$contrasts4.3.1.9 Visualize the model estimates
plot(conditional_effects(mnb_iconic_cond_round_log), ask = FALSE)
5 =====<Chosen> Rate of iconic gestures per 100 words=====
5.1 —DataViz: rate per 100 words—
5.1.1 bp: mean by condition
bp_mean_iconic_rate_by_cond = ggplot(data=trial_info_pair,
aes(x=condition, y=mean_iconic_per_100words, fill=condition)) +
geom_jitter(aes(color = pair),
size = 1, alpha = 1,
width = 0.07, height = 0) +
geom_boxplot(width = .2,
outlier.shape = NA, alpha = 0.7) +
geom_point(data = trial_info_cond,
aes(y = mean_iconic_per_100words),
shape = 21, size = 3, fill = "white") +
scale_y_continuous(limits = c(0, 11), breaks = seq(0, 10, 2)) +
scale_fill_manual(values = c("#ED6B06", "#00786A", "darkgrey")) +
labs(x = "Visibility",
y = "Mean rate of iconic gestures") +
theme_clean(base_size = 15) +
theme(axis.text.x = element_text(colour = "black", size = 14),
axis.text.y = element_text(colour = "black", size = 14),
axis.title = element_text(size = 15, face = 'bold'),
axis.title.x = element_text(vjust = -2),
axis.title.y = element_text(vjust = 2),
legend.position = "none",
strip.text = element_text(size = 15, face = 'bold'),
plot.background = element_blank(),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines"))
ggplotly(bp_mean_iconic_rate_by_cond)5.1.2 bp: mean by condition x round
pd = position_dodge(width = .75)
ggplot(data=df_trial_info,
aes(x=round, y=n_iconic_per_100words, fill=condition)) +
geom_boxplot(outlier.shape = NA,
alpha = 0.7) +
geom_point(data = trial_info_cond_round,
aes(y = mean_iconic_per_100words,
group = condition),
position = pd,
shape = 21, size = 2, fill = "white") +
scale_fill_manual(values = c("#ED6B06", "#00786A", "darkgrey")) +
scale_y_continuous(limits = c(0, 22), breaks = seq(0, 20, 5)) +
labs(x = "Round",
y = "Mean rate of iconic gestures") +
theme_clean(base_size = 15) +
theme(axis.text.x = element_text(colour = "black", size = 14),
axis.text.y = element_text(colour = "black", size = 14),
axis.title = element_text(size = 15, face = 'bold'),
axis.title.x = element_text(vjust = -2),
axis.title.y = element_text(vjust = 2),
legend.position = "top",
strip.text = element_text(size = 15, face = 'bold'),
plot.background = element_blank(),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines"))
5.2 —<Chosen> ZI negative binomial regression models—
5.2.1 Prior specification
As the unit of the dependent variable is different from the previous model, we will set different priors for the rate of iconic gestures per 100 words. We will use weakly informative priors for the regression each parameter. To specify priors for the intercept, we will refer to the rate of iconic gestures per 100 words reported in Akamine et al. (2024). In the paper, the authors analyzed data from 19 dyads involving in the same task as the current study but in co-present interaction. The total number of iconic gestures reported was 4413, which was collected from 19 dyads, each performing 96 trials. The total number of words produced was 71695 (28152 content) words. Therefore, the expected rate of iconic gestures per 100 words is \(4413 / (71695 / 100) = 6.16\) (log(6.16) = 1.82).
For the fixed effects, we will set unbiased priors with a mean of 0 and a standard deviation of 0.5.
For the standard deviation of the random effects, we set the prior to be normal with mean 0 and standard deviation 0.5. For the correlation between the random effects, we set the prior to be LKJ(2).
For the models including the round as fixed effects (whether backward-difference coded or centered + log-transformed), the intercept will represent the mean expected number of iconic gestures (ground mean). As the meaning of the intercept doesn’t change when adding the round variable, we use the same prior for the intercept.
Note that we do not expect the rate of iconic gestures to change across rounds (i.e., we expect the number of words and gestures to decrease at an approximately same pace across the rounds). Also, it is common to set the mean of slopes to be 0 to avoid favoring any directions. Therefore we will set the mean of the prior for slope to 0.
priors_rslope_rate = c(
prior(normal(1.82, 0.5), class = Intercept),
prior(normal(0, 0.5), class = b),
prior(normal(0, 0.5), class = sd),
prior(lkj(2), class = cor))
priors_rslope_rate_zinb = c(
prior(normal(1.82, 0.5), class = Intercept),
prior(normal(0, 0.5), class = b),
prior(normal(0, 0.5), class = sd),
prior(lkj(2), class = cor),
prior(normal(0, 0.5), class = zi), # on the logit scale
prior(normal(0, 50), class = shape))5.2.2 Model 3: condition x round (pair)
In the previous section, we analyzed the number of iconic gestures produced per trial. However, it is common to analyze the rate of iconic gestures per 100 words to account for the differences in the length of the trials and speech rate. Here, we will include the log of speech rate (number of words / 100) as an exposure variable and analyze the rate of iconic gestures per 100 words by condition. Note that the syntax for the exposure variable is different from the Poisson regression model; for negative binomial regression, the exposure variable is included with the rate() function.
5.2.2.1 Model Effect: round
nb_iconic_rate_cond_round = brm(num_iconic_gestures | rate(n_words_100) ~
1 + condition * round +
(1+round|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope_rate,
data = df_trial_info,
sample_prior = T,
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/nb_iconic_rate_cond_round")
nb_iconic_rate_cond_round_c = brm(num_iconic_gestures | rate(n_words_100) ~
1 + condition * round_c +
(1+round_c|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope_rate,
data = df_trial_info,
sample_prior = T,
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/nb_iconic_rate_cond_round_c")
nb_iconic_rate_cond_round_log = brm(num_iconic_gestures | rate(n_words_100) ~
1 + condition * log_round_c +
(1+log_round_c|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope_rate,
data = df_trial_info,
sample_prior = T,
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/nb_iconic_rate_cond_round_log")
### loo compare
if (!file.exists("models/loo_nb_iconic_rate_cond_round.rds")){
nb_cond_round_loo = loo(nb_iconic_rate_cond_round)
saveRDS(nb_cond_round_loo, file = "models/loo_nb_iconic_rate_cond_round.rds")
}
if (!file.exists("models/loo_nb_iconic_rate_cond_round_c.rds")){
nb_cond_round_c_loo = loo(nb_iconic_rate_cond_round_c)
saveRDS(nb_cond_round_c_loo, file = "models/loo_nb_iconic_rate_cond_round_c.rds")
}
if (!file.exists("models/loo_nb_iconic_rate_cond_round_log.rds")){
nb_cond_round_log_loo = loo(nb_iconic_rate_cond_round_log)
saveRDS(nb_cond_round_log_loo, file = "models/loo_nb_iconic_rate_cond_round_log.rds")
}
nb_cond_round_loo = readRDS("models/loo_nb_iconic_rate_cond_round.rds")
nb_cond_round_c_loo = readRDS("models/loo_nb_iconic_rate_cond_round_c.rds")
nb_cond_round_log_loo = readRDS("models/loo_nb_iconic_rate_cond_round_log.rds")
loo_compare(nb_cond_round_loo, nb_cond_round_c_loo, nb_cond_round_log_loo)## elpd_diff se_diff
## nb_iconic_rate_cond_round_c 0.0 0.0
## nb_iconic_rate_cond_round -9.8 5.3
## nb_iconic_rate_cond_round_log -13.0 4.2The leave-one-out (LOO) Effect shows that centered round provides a larger predictive power (elpd_diff) and smaller standard error (se_diff) compared to the backward-difference coded round or centered log-transformed round. Therefore, we will use the centered round.
5.2.2.2 Model Effect: ZI or not
zinb_iconic_rate_cond_round_c = brm(num_iconic_gestures ~ 1 + condition * round_c +
offset(log(n_words_100)) +
(1+round_c|pair) + (1|target),
family = zero_inflated_negbinomial(),
prior = priors_rslope_rate_zinb,
data = df_trial_info,
sample_prior = T,
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/zinb_iconic_rate_cond_round_c")
### loo compare
if (!file.exists("models/loo_zinb_iconic_rate_cond_round_c.rds")){
zinb_cond_round_c_loo = loo(zinb_iconic_rate_cond_round_c)
saveRDS(zinb_cond_round_c_loo, file = "models/loo_zinb_iconic_rate_cond_round_c.rds")
}
nb_cond_round_c_loo = readRDS("models/loo_nb_iconic_rate_cond_round_c.rds")
zinb_cond_round_c_loo = readRDS("models/loo_zinb_iconic_rate_cond_round_c.rds")
loo_compare(nb_cond_round_c_loo, zinb_cond_round_c_loo)## elpd_diff se_diff
## zinb_iconic_rate_cond_round_c 0.0 0.0
## nb_iconic_rate_cond_round_c -13.7 6.2The leave-one-out (LOO) Effect shows that zero-inflation model has a higher predictive power. As such, we will use zero-inflated negative binomial regression model.
5.2.2.3 Prior predictive check
zinb_iconic_rate_cond_round_c_prior = brm(num_iconic_gestures ~ 1 + condition * round_c +
offset(log(n_words_100)) +
(1+round_c|pair) + (1|target),
family = zero_inflated_negbinomial(),
prior = priors_rslope_rate_zinb,
sample_prior = "only",
data = df_trial_info,
file = "models/zinb_iconic_rate_cond_round_c_prior")
pp_check(zinb_iconic_rate_cond_round_c_prior, ndraws = 100, type = "bars") +
coord_cartesian(xlim = c(0, 20),
ylim = c(0, 4000))
The prior predictive check shows that the model generates data that are somewhat similar to the observed data.
5.2.2.4 Fit the model
zinb_iconic_rate_cond_round_c = brm(num_iconic_gestures ~ 1 + condition * round_c +
offset(log(n_words_100)) +
(1+round_c|pair) + (1|target),
family = zero_inflated_negbinomial(),
prior = priors_rslope_rate_zinb,
data = df_trial_info,
sample_prior = T,
save_pars = save_pars(all = TRUE),
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/zinb_iconic_rate_cond_round_c")
model = zinb_iconic_rate_cond_round_c
summary(model)## Family: zero_inflated_negbinomial
## Links: mu = log; shape = identity; zi = identity
## Formula: num_iconic_gestures ~ 1 + condition * round_c + offset(log(n_words_100)) + (1 + round_c | pair) + (1 | target)
## Data: df_trial_info (Number of observations: 4316)
## Draws: 4 chains, each with iter = 20000; warmup = 2000; thin = 1;
## total post-warmup draws = 72000
##
## Multilevel Hyperparameters:
## ~pair (Number of levels: 45)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.66 0.07 0.53 0.82 1.00 18554
## sd(round_c) 0.11 0.02 0.08 0.15 1.00 26924
## cor(Intercept,round_c) 0.71 0.11 0.46 0.87 1.00 38762
## Tail_ESS
## sd(Intercept) 34315
## sd(round_c) 44094
## cor(Intercept,round_c) 50332
##
## ~target (Number of levels: 16)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.10 0.03 0.06 0.16 1.00 26905 44443
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 1.33 0.10 1.13 1.54 1.00 12630
## conditionAO_Asym 0.28 0.22 -0.15 0.70 1.00 14119
## conditionAsym_Sym -0.06 0.22 -0.49 0.36 1.00 13682
## round_c -0.13 0.02 -0.17 -0.09 1.00 21391
## conditionAO_Asym:round_c 0.02 0.05 -0.07 0.11 1.00 24222
## conditionAsym_Sym:round_c -0.03 0.04 -0.12 0.06 1.00 22801
## Tail_ESS
## Intercept 23347
## conditionAO_Asym 25910
## conditionAsym_Sym 24753
## round_c 36599
## conditionAO_Asym:round_c 40686
## conditionAsym_Sym:round_c 36775
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 41.96 16.08 21.47 83.70 1.00 100817 52037
## zi 0.04 0.01 0.02 0.06 1.00 99853 44866
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).bayestestR::hdi(model)# bayestestR::hdi(model, ci = 0.89)The coefficients show that the condition does not have a significant effect on the rate of iconic gestures per 100 words. However, there was a significant decrease in the rate of iconic gestures per 100 words across the rounds. This means that the number of iconic gestures decreased more than the number of words did across the rounds. A formal hypothesis testing will be performed later using Bayes factor.
5.2.2.5 Visualize the posterior distributions
# p = plot_posterior(model, interaction = T)
# p5.2.2.6 Hypothesis testing: Bayes factor
### varying priors for sensitivity analysis
prior_size = c("xs", "s", "l", "xl")
prior_sd = c(0.1, 0.3, 1, 1.5)
bfs_cond_ao_asym = c()
bfs_cond_asym_sym = c()
bfs_round = c()
for (i in 1:length(prior_sd)){
priors = c(
prior(normal(1.82, 0.5), class = Intercept),
set_prior(paste0("normal(0,", prior_sd[i], ")"), class = "b"),
prior(normal(0, 0.5), class = sd),
prior(lkj(2), class = cor),
prior(normal(0, 0.5), class = zi), # on the logit scale
prior(normal(0, 50), class = shape))
fname = paste0("models/zinb_iconic_rate_cond_round_c_", prior_size[i])
fit = brm(num_iconic_gestures ~
1 + condition * round_c +
offset(log(n_words_100)) +
(1+round_c|pair) + (1|target),
family = zero_inflated_negbinomial(),
prior = priors,
data = df_trial_info,
sample_prior = T,
save_pars = save_pars(all = TRUE),
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = fname)
# BF for sym - asym
h = hypothesis(fit, "conditionAO_Asym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_cond_ao_asym = c(bfs_cond_ao_asym, bf)
# BF for sym - ao
h = hypothesis(fit, "conditionAsym_Sym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_cond_asym_sym = c(bfs_cond_asym_sym, bf)
# BF for round
h = hypothesis(fit, "round_c = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round = c(bfs_round, bf)
}
### add BF for the main/medium model
prior_size[3:5] = c("m", prior_size[3:4])
prior_sd[3:5] = c(0.5, prior_sd[3:4])
# BF for sym - asym
h = hypothesis(model, "conditionAO_Asym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_cond_ao_asym[3:5] = c(bf, bfs_cond_ao_asym[3:4])
# BF for sym - ao
h = hypothesis(model, "conditionAsym_Sym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_cond_asym_sym[3:5] = c(bf, bfs_cond_asym_sym[3:4])
# BF for round
h = hypothesis(model, "round_c = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round[3:5] = c(bf, bfs_round[3:4])
### make a df for BFs
df_bf = data.frame(size = prior_size,
sd = prior_sd,
ao_asym = bfs_cond_ao_asym,
asym_sym = bfs_cond_asym_sym,
round = bfs_round) %>%
mutate(prior = paste0("N(0, ", sd, ")")) %>%
pivot_longer(cols = c("ao_asym", "asym_sym", "round"),
names_to = "Effect",
values_to = "BF10") %>%
mutate(Effect = factor(Effect,
levels = c("ao_asym", "asym_sym", "round"),
labels = c("AO-AsymAV", "AsymAV-SymAV", "Round")),
Predictor = ifelse(Effect == "round", "Round", "Visibility"),
BF10_log10 = log10(BF10))
df_bf %>% arrange(Effect, sd)#### Plot BFs ####
ggplot(df_bf, aes(x = prior, y = BF10, group = Effect)) +
geom_hline(yintercept = 1, linetype="dashed") +
geom_point(aes(color=Effect)) +
geom_line(aes(color=Effect)) +
facet_wrap(vars(Predictor), scales="free_y") +
theme_bw(base_size = 12)+
theme(legend.position = "top")+
scale_y_log10("Bayes factor (BF10)",
breaks = c(0.001, 0.03, 0.01, 0.1, 0.33, 1, 3, 10, 30, 100),
labels = c(0.001, 0.03, 0.01, 0.1, 0.33, 1, 3, 10, 30, 100)) +
xlab("prior")
5.2.2.7 Posterior predictive check
model = zinb_iconic_rate_cond_round_c
pp_check_sym = pp_check_each_condition(model, df_trial_info, "SymAV")
pp_check_asym = pp_check_each_condition(model, df_trial_info, "AsymAV")
pp_check_ao = pp_check_each_condition(model, df_trial_info, "AO")
gridExtra::grid.arrange(pp_check_sym, pp_check_asym, pp_check_ao, ncol = 2)
pp_check_r1 = pp_check_each_round(model, df_trial_info, "R1")
pp_check_r2 = pp_check_each_round(model, df_trial_info, "R2")
pp_check_r3 = pp_check_each_round(model, df_trial_info, "R3")
pp_check_r4 = pp_check_each_round(model, df_trial_info, "R4")
pp_check_r5 = pp_check_each_round(model, df_trial_info, "R5")
pp_check_r6 = pp_check_each_round(model, df_trial_info, "R6")
gridExtra::grid.arrange(pp_check_r1, pp_check_r2, pp_check_r3, pp_check_r4, pp_check_r5, pp_check_r6, ncol = 3)
5.2.2.8 Pair-wise Effect
emmeans(model, pairwise ~ condition)$contrasts## contrast estimate lower.HPD upper.HPD
## SymAV - AsymAV -0.0633 -0.497 0.355
## SymAV - AO 0.2134 -0.242 0.678
## AsymAV - AO 0.2760 -0.140 0.708
##
## Point estimate displayed: median
## Results are given on the log (not the response) scale.
## HPD interval probability: 0.95emmeans(model, pairwise ~ condition, level = 0.89)$contrasts## contrast estimate lower.HPD upper.HPD
## SymAV - AsymAV -0.0633 -0.405 0.286
## SymAV - AO 0.2134 -0.163 0.584
## AsymAV - AO 0.2760 -0.072 0.616
##
## Point estimate displayed: median
## Results are given on the log (not the response) scale.
## HPD interval probability: 0.895.2.2.9 Visualize the model estimates
plot(conditional_effects(model), ask = FALSE)
6 =====Causal model of gestural alignment=====
Some experts in the field of statistics and causal inference have advised that researchers should build a causal model and examine which factors should be included and excluded from regression models if their aim is to infer the causal effects (e.g., McElreath, 2020, Pearl, 2010, Cinelli, Forney, & Pearl, 2020). Following this advice, we will build a causal model to infer the direct effect of speaker visibility on gestural alignment.
We assume the following causal model:
### causal model for gestural alignment rate
dag_gest <- dagitty('dag {
visibility -> gest_align
visibility -> n_words
visibility -> n_gestures
n_words -> lex_align
n_words -> n_gestures
round -> n_words
round -> lex_align
round -> gest_align
round -> n_gestures
n_gestures -> gest_align
lex_align -> gest_align
}')
plot(dag_gest)
### check the minimam adjustment set
print("Direct effect of visibility on gestural alignment rate")## [1] "Direct effect of visibility on gestural alignment rate"adjustmentSets(dag_gest, exposure = "visibility", outcome = "gest_align", effect = "direct")## { lex_align, n_gestures, round }
## { n_gestures, n_words, round }print("Direct effect of lexical alignment rate on gestural alignment rate")## [1] "Direct effect of lexical alignment rate on gestural alignment rate"adjustmentSets(dag_gest, exposure = "lex_align", outcome = "gest_align", effect = "direct")## { n_gestures, round, visibility }
## { n_words, round }print("Direct effect of round on gestural alignment rate")## [1] "Direct effect of round on gestural alignment rate"adjustmentSets(dag_gest, exposure = "round", outcome = "gest_align", effect = "direct")## { lex_align, n_gestures, visibility }The minimum adjustment set for estimating the direct causal effect of speaker visibility on gestural alignment rate is {lex_align, n_words, round}. Note that because dagitty didn’t find any adjustment set for the direct effect of visibility on lexical alignment rate when we expected bidirectional causation between lexical and gestural alignment, we assumed a unidirecitonal causation from lexical alignment to gestural alignment only in this model. This will be reversed in the causal model for lexical alignment rate, such that we assume a unidirectional causation from gestural alignment to lexical alignment.
In addition, we are also interested in the direct effect of lexical alignment rate on gestural alignment rate. The minimum adjustment set for is {visibility, n_gestures, round}.
As the minimum adjustment sets for the direct effects of visibility, lexical alignment rate, and round on gestural alignment rate are identical (i.e., {visibility, lex_align, n_gestures, round}), we can estimate the direct effect of these variables on gestural alignment rate with the following model:
\[ gest\_align \: | \: \text{rate}(n\_iconic\_gesture) \sim \\ \text{visibility} * \text{round} + \text{n_lexical_alignment} + \\ (1+\text{round} | \text{participant}) + (1 | \text{item}) \]
Just to be thorough, we will first analyze the raw frequency and rate
of gestural alignment. If you just want to see the analysis of the
proportion of gestural alignment, see [-–
7 =====Number of gestural alignment=====
7.1 —DataViz: number of gestural alignment—
7.1.1 bp: mean by condition
bp_mean_gest_alignment_by_cond = ggplot(data=trial_info_pair,
aes(x=condition, y=mean_num_gest_align, fill=condition)) +
geom_jitter(aes(color = pair),
size = 1, alpha = 1,
width = 0.07, height = 0) +
geom_boxplot(width = .4,
outlier.shape = NA, alpha = 0.7) +
geom_point(data = trial_info_cond,
aes(y = mean_num_gest_align),
size = 2, shape = 4) +
scale_fill_manual(values = c("#ED6B06", "#00786A", "darkgrey")) +
labs(x = "Visibility",
y = "Mean gestural alignment count") +
theme_clean(base_size = 15) +
theme(axis.text.x = element_text(colour = "black", size = 14),
axis.text.y = element_text(colour = "black", size = 14),
axis.title = element_text(size = 15, face = 'bold'),
axis.title.x = element_text(vjust = -2),
axis.title.y = element_text(vjust = 2),
legend.position = "none",
strip.text = element_text(size = 15, face = 'bold'),
plot.background = element_blank(),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines"))
ggplotly(bp_mean_gest_alignment_by_cond)7.2 —Negative binomial regression—
7.2.1 Prior specification
We will set priors for the intercept based on the expected number of gestural alignment reported in Akamine et al. (2024). In the paper, the authors analyzed data from 19 dyads involving in the same task as the current study but in co-present interaction. The total number of gestural alignment reported was 1086, which was collected from 19 dyads, each performing 96 trials. This leads to the expected number of gestural alignment is 1086 / (19 * 96) = 0.6 (log(0.6) = -0.5). As such, we will set the prior for the intercept to be normal with a mean of -0.5 and a standard deviation of 0.5.
For the fixed effects, we will set unbiased priors with a mean of 0. We set different SDs for each effect because the scale of each predictor is different. For N. lexical alignment and N. iconic gestures, we set a prior of Normal(0, 0.2). This means that if the mean N. lex align or N. iconic gesture increase by 1, we expect the mean N. gestural alignment to change by -0.6 to 0.9, with the most likely size of change to be 0. As for the visibility condition and the round, we set a prior of Normal(0. 0.5).
- N. lex align: Normal(0, 0.2) <– range of expected effect size in the original scale is [-0.6, 0.9]
- N. iconic gestures: Normal(0, 0.2) <– [-0.6, 0.9]
- visibility condition: Normal(0, 0.5) <– [-1.27, 3.13]
- Round: Normal(0, 0.5) <– [-1.27, 3.13]
We will also specify a prior for the shape parameter to prevent the model from returning divergent transitions. We will set the prior to be normal with a mean of 0 and a wide standard deviation of 50.
For the standard deviation of the random effects, we set the prior to be normal with mean 0 and standard deviation 0.5. For the correlation between the random effects, we set the prior to be LKJ(2).
### pair
priors_rint_gest_align = c(
prior(normal(0.6, 0.5), class = Intercept),
prior(normal(0, 0.5), class = b),
prior(normal(0, 0.2), class = b, coef = "num_lex_align"),
prior(normal(0, 0.2), class = b, coef = "num_iconic_gestures"),
prior(normal(0, 50), class = shape),
prior(normal(0, 0.5), class = sd))
priors_rslope_gest_align = c(
prior(normal(0.6, 0.5), class = Intercept),
prior(normal(0, 0.5), class = b),
prior(normal(0, 0.2), class = b, coef = "num_lex_align"),
prior(normal(0, 0.2), class = b, coef = "num_iconic_gestures"),
prior(normal(0, 50), class = shape),
prior(normal(0, 0.5), class = sd),
prior(lkj(2), class = cor))
priors_rslope_gest_align_zinb = c(
prior(normal(0.6, 0.5), class = Intercept),
prior(normal(0, 0.5), class = b),
prior(normal(0, 0.2), class = b, coef = "num_lex_align"),
prior(normal(0, 0.2), class = b, coef = "num_iconic_gestures"),
prior(normal(0, 0.5), class = sd),
prior(lkj(2), class = cor),
prior(normal(0, 0.5), class = zi), # on the logit scale
prior(normal(0, 50), class = shape))7.2.2 Model 4: [pair] effect of lex align on gest align
7.2.2.1 Prior predictive check
zinb_gest_align_prior = brm(num_gestural_align ~
1 + condition + round + num_lex_align + num_iconic_gestures +
(1+round|pair) + (1|target),
family = zero_inflated_negbinomial(),
prior = priors_rslope_gest_align_zinb,
data = df_trial_info,
sample_prior = "only",
control = list(adapt_delta = 0.9,
max_treedepth = 20),
file = "models/zinb_gest_align_prior")
pp_check(zinb_gest_align_prior, ndraws = 100, type = "bars") +
coord_cartesian(xlim = c(0, 20),
ylim = c(0, 4000))
The prior predictive check shows that the model expects fewer amount of 2, 3, and 4. This suggests that the zero-inflation prior may be too large or the mean for the intercept prior is too low. We will check the prior-posterior update plot and posterior predictive check to see if the model generates data that are similar to the observed data. If not, we will consider modifying the priors.
7.2.2.2 Fit the model
zinb_align_cond_round = brm(num_gestural_align ~
1 + condition + round + num_lex_align + num_iconic_gestures +
(1+round|pair) + (1|target),
family = zero_inflated_negbinomial(),
prior = priors_rslope_gest_align_zinb,
data = df_trial_info,
sample_prior = T,
save_pars = save_pars(all = TRUE),
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/zinb_align_cond_round")
model = zinb_align_cond_round
summary(model)## Family: zero_inflated_negbinomial
## Links: mu = log; shape = identity; zi = identity
## Formula: num_gestural_align ~ 1 + condition + round + num_lex_align + num_iconic_gestures + (1 + round | pair) + (1 | target)
## Data: df_trial_info (Number of observations: 4316)
## Draws: 4 chains, each with iter = 20000; warmup = 2000; thin = 1;
## total post-warmup draws = 72000
##
## Multilevel Hyperparameters:
## ~pair (Number of levels: 45)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 1.04 0.12 0.82 1.30 1.00 13258
## sd(roundR12) 0.62 0.15 0.33 0.93 1.00 24974
## sd(roundR23) 0.29 0.12 0.05 0.52 1.00 14676
## sd(roundR34) 0.19 0.12 0.01 0.46 1.00 12288
## sd(roundR45) 0.22 0.13 0.01 0.49 1.00 16214
## sd(roundR56) 0.16 0.12 0.01 0.45 1.00 28127
## cor(Intercept,roundR12) 0.33 0.21 -0.12 0.69 1.00 38507
## cor(Intercept,roundR23) 0.47 0.25 -0.13 0.85 1.00 36116
## cor(roundR12,roundR23) 0.06 0.27 -0.47 0.59 1.00 43615
## cor(Intercept,roundR34) 0.18 0.31 -0.48 0.72 1.00 44648
## cor(roundR12,roundR34) 0.10 0.31 -0.52 0.66 1.00 47546
## cor(roundR23,roundR34) 0.16 0.31 -0.48 0.72 1.00 41117
## cor(Intercept,roundR45) 0.15 0.31 -0.48 0.69 1.00 50329
## cor(roundR12,roundR45) 0.16 0.30 -0.46 0.70 1.00 45069
## cor(roundR23,roundR45) 0.20 0.31 -0.46 0.74 1.00 31625
## cor(roundR34,roundR45) 0.11 0.32 -0.54 0.69 1.00 35440
## cor(Intercept,roundR56) -0.08 0.33 -0.67 0.57 1.00 68954
## cor(roundR12,roundR56) 0.04 0.32 -0.58 0.64 1.00 62152
## cor(roundR23,roundR56) 0.04 0.33 -0.60 0.65 1.00 60123
## cor(roundR34,roundR56) 0.05 0.33 -0.59 0.66 1.00 50904
## cor(roundR45,roundR56) 0.03 0.33 -0.60 0.65 1.00 53360
## Tail_ESS
## sd(Intercept) 27474
## sd(roundR12) 27105
## sd(roundR23) 11704
## sd(roundR34) 20523
## sd(roundR45) 19835
## sd(roundR56) 30531
## cor(Intercept,roundR12) 46914
## cor(Intercept,roundR23) 32737
## cor(roundR12,roundR23) 47771
## cor(Intercept,roundR34) 48873
## cor(roundR12,roundR34) 53571
## cor(roundR23,roundR34) 49423
## cor(Intercept,roundR45) 48204
## cor(roundR12,roundR45) 49810
## cor(roundR23,roundR45) 44304
## cor(roundR34,roundR45) 48265
## cor(Intercept,roundR56) 53547
## cor(roundR12,roundR56) 53118
## cor(roundR23,roundR56) 57059
## cor(roundR34,roundR56) 56039
## cor(roundR45,roundR56) 56078
##
## ~target (Number of levels: 16)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.22 0.06 0.13 0.35 1.00 21997 37924
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -1.93 0.17 -2.26 -1.58 1.00 7722 15921
## conditionAO_Asym 0.32 0.29 -0.25 0.88 1.00 13178 24373
## conditionAsym_Sym 0.25 0.27 -0.29 0.79 1.00 11227 21534
## roundR12 1.54 0.16 1.23 1.85 1.00 31455 40797
## roundR23 -0.10 0.11 -0.32 0.10 1.00 30196 40085
## roundR34 -0.25 0.11 -0.47 -0.06 1.00 35515 38965
## roundR45 -0.31 0.12 -0.56 -0.08 1.00 39319 44558
## roundR56 -0.18 0.12 -0.43 0.05 1.00 58580 50167
## num_lex_align 0.05 0.02 0.00 0.09 1.00 68008 56486
## num_iconic_gestures 0.20 0.01 0.18 0.23 1.00 41732 49580
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 3.65 0.53 2.77 4.83 1.00 49170 51339
## zi 0.01 0.01 0.00 0.03 1.00 72219 39837
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).bayestestR::hdi(model)The coefficients show that the number of lexical alignment is a significant predictor of the number of gestural alignment.
7.2.2.3 Hypothesis testing: Bayes factor
### varying priors for sensitivity analysis
prior_size = c("xs", "s", "l", "xl")
prior_sd = c(0.05, 0.1, 0.3, 0.5)
bfs_lex_align = c()
for (i in 1:length(prior_sd)){
priors = c(
prior(normal(0.6, 0.5), class = Intercept),
prior(normal(0, 0.5), class = b),
set_prior(paste0("normal(0,", prior_sd[i], ")"), class = "b", coef = "num_lex_align"),
prior(normal(0, 0.2), class = b, coef = "num_iconic_gestures"),
prior(normal(0, 0.5), class = sd),
prior(lkj(2), class = cor),
prior(normal(0, 0.5), class = zi),
prior(normal(0, 50), class = shape))
fname = paste0("models/zinb_align_cond_round_", prior_size[i])
fit = brm(num_gestural_align ~
1 + condition + round + num_lex_align + num_iconic_gestures +
(1+round|pair) + (1|target),
family = zero_inflated_negbinomial(),
prior = priors,
data = df_trial_info,
sample_prior = T,
save_pars = save_pars(all = TRUE),
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = fname)
### BF for N. lex alignment
h = hypothesis(fit, "num_lex_align = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_lex_align = c(bfs_lex_align, bf)
}
### add BF for the main/medium model
prior_size[3:5] = c("m", prior_size[3:4])
prior_sd[3:5] = c(0.2, prior_sd[3:4])
### BF for N. lex alignment
h = hypothesis(model, "num_lex_align = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_lex_align[3:5] = c(bf, bfs_lex_align[3:5])
### make a df for BFs
df_bf_lex = data.frame(size = prior_size,
sd = prior_sd,
lex_align = bfs_lex_align) %>%
mutate(prior = paste0("N(0, ", sd, ")")) %>%
pivot_longer(cols = c("lex_align"),
names_to = "Effect",
values_to = "BF10") %>%
mutate(Effect = "N. lex align",
Predictor = "N. lex align")#### Plot BFs ####
ggplot(filter(df_bf_lex),
aes(x = factor(sd), y = BF10, group = Effect)) +
geom_hline(yintercept = 1, linetype="dashed") +
geom_point(aes(color=Effect)) +
geom_line(aes(color=Effect)) +
facet_wrap(vars(Predictor)) +
theme_clean(base_size = 15) +
theme(axis.text.x = element_text(colour = "black", size = 14),
axis.text.y = element_text(colour = "black", size = 14),
axis.title = element_text(size = 15, face = 'bold'),
axis.title.x = element_text(vjust = -2),
axis.title.y = element_text(vjust = 2),
legend.position = "top",
strip.text = element_text(size = 15, face = 'bold'),
plot.background = element_blank(),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines")) +
scale_y_log10("Bayes factor (BF10)",
breaks = c(0.001, 0.03, 0.01, 0.1, 0.33, 1, 3, 10, 30, 100),
labels = c(0.001, 0.03, 0.01, 0.1, 0.33, 1, 3, 10, 30, 100)) +
xlab("SD for the prior")
8 =====<Chosen> Proportion of gestural alignment (n_gest_alignment / n_iconic)=====
We can analyze the proportion of gestural alignment in two ways: (1) modeling the rate of gestural alignment per iconic gesture using a negative binomial regression model and (2) modeling the probability of gestural alignment using a zero-one-inflated beta regression model.
Key differences in the two models are that the negative binomial regression model assumes that the rate of gestural alignment is a count variable, while the zero-one-inflated beta regression model assumes that the proportion of gestural alignment is a continuous variable bounded between 0 and 1. Also, while negative binomial regression models the rate of events, the zero-one-inflated beta regression models the probability of events. In the case of the proportion of gestural alignment, two models should yield similar results, but it is important to note that the two models are conceptually different.
8.1 —DataViz: proportion of gestural alignment—
8.1.1 bp: mean by condition
bp_mean_gest_alignment_prop_by_cond = ggplot(data=trial_info_pair,
aes(x=condition, y=mean_gest_align_prop, fill=condition)) +
geom_jitter(aes(color = pair),
size = 1, alpha = 1,
width = 0.07, height = 0) +
geom_boxplot(width = .3,
outlier.shape = NA, alpha = 0.7) +
geom_point(data = trial_info_cond,
aes(y = mean_gest_align_prop),
shape = 21, size = 3, fill = "white") +
scale_fill_manual(values = c("#ED6B06", "#00786A", "darkgrey")) +
labs(x = "Visibility",
y = "Mean gest align rate") +
theme_clean(base_size = 15) +
theme(axis.text.x = element_text(colour = "black", size = 14),
axis.text.y = element_text(colour = "black", size = 14),
axis.title = element_text(size = 15, face = 'bold'),
axis.title.x = element_text(vjust = -2),
axis.title.y = element_text(vjust = 2),
legend.position = "none",
strip.text = element_text(size = 15, face = 'bold'),
plot.background = element_blank(),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines"))
ggplotly(bp_mean_gest_alignment_prop_by_cond)8.1.2 bp: mean by condition and round
pd = position_dodge(width = .75)
bp_mean_gest_alignment_prop_by_round_cond =
ggplot(data=trial_info_pair_round,
aes(x=round_n, y=mean_gest_align_prop, fill = condition)) +
geom_jitter(aes(color = pair),
size = 0.5, alpha = 0.7,
width = 0.07, height = 0) +
geom_boxplot(width = .5,
outlier.shape = NA, alpha = 0.7) +
geom_point(data = trial_info_cond_round,
aes(y = mean_gest_align_prop),
position = pd,
shape = 21, size = 2, fill = "white") +
scale_fill_manual(values = c("#ED6B06", "#00786A", "darkgrey")) +
scale_y_continuous(limits = c(0, 1)) +
labs(x = "Round",
y = "Mean gest align rate") +
theme_clean(base_size = 15) +
theme(axis.text.x = element_text(colour = "black", size = 14),
axis.text.y = element_text(colour = "black", size = 14),
axis.title = element_text(size = 15, face = 'bold'),
axis.title.x = element_text(vjust = -2),
axis.title.y = element_text(vjust = 2),
legend.position = "none",
strip.text = element_text(size = 15, face = 'bold'),
plot.background = element_blank(),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines")) +
facet_grid(cols = vars(condition))
ggplotly(bp_mean_gest_alignment_prop_by_round_cond)8.2 Prepare data
To model the proportion of gestural alignment, we need to remove trials where the number of iconic gestures is 0. This is because dividing any numbers by 0 results in an undefined value (NA) and prevents model from running.
df_gest_align_posreg_prop = df_trial_info %>%
filter(num_iconic_gestures > 0)
print(paste0("Number of rows before removing trials with no iconic gestures: ", nrow(df_trial_info)))## [1] "Number of rows before removing trials with no iconic gestures: 4316"print(paste0("Number of rows before after trials with no iconic gestures: ", nrow(df_gest_align_posreg_prop)))## [1] "Number of rows before after trials with no iconic gestures: 2199"print(paste0("Number of removed trials: ", nrow(df_trial_info) - nrow(df_gest_align_posreg_prop)))## [1] "Number of removed trials: 2117"8.3 Prior specification
We will set priors based on Akamine et al. (2024). As mentioned in the previous analysis, they detected 4413 iconic gestures and 1086 instances of gestural alignment. Dividing the number of gestural alignment by the number of iconic gestures gives 0.25 (1086/4413) (-1.39 on the log scale). This means that we expect 1 gestural alignment per 4 iconic gestures.
For the slope parameters, we set the mean of 0 with a SD of 0.5. This means that for example, we expect the mean difference between the AO and AsymAV conditions to range from -0.16 to 0.43.
### pair
priors_rint_gest_align_prop = c(
prior(normal(-1.39, 0.5), class = Intercept),
prior(normal(0, 0.5), class = b),
prior(normal(0, 50), class = shape),
prior(normal(0, 0.5), class = sd))
priors_rslope_gest_align_prop = c(
prior(normal(-1.39, 0.5), class = Intercept),
prior(normal(0, 0.5), class = b),
prior(normal(0, 50), class = shape),
prior(normal(0, 0.5), class = sd),
prior(lkj(2), class = cor))
priors_rslope_gest_align_prop_zinb = c(
prior(normal(-1.39, 0.5), class = Intercept),
prior(normal(0, 0.5), class = b),
prior(normal(0, 0.5), class = sd),
prior(lkj(2), class = cor),
prior(normal(0, 0.5), class = zi), # on the logit scale
prior(normal(0, 50), class = shape))8.4 —Negative binomial regression models—
8.4.1 Model 5: [pair] condition * round + num_lex_align
8.4.1.1 Model Effect: round
nb_align_rate_cond_round = brm(num_gestural_align | rate(num_iconic_gestures) ~
1 + condition * round + num_lex_align +
(1+round|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope_gest_align_prop,
data = df_gest_align_posreg_prop,
sample_prior = T,
save_pars = save_pars(all = TRUE),
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/nb_align_rate_cond_round")
nb_align_rate_cond_round_c = brm(num_gestural_align | rate(num_iconic_gestures) ~
1 + condition * round_c + num_lex_align +
(1+round_c|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope_gest_align_prop,
data = df_gest_align_posreg_prop,
sample_prior = T,
save_pars = save_pars(all = TRUE),
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/nb_align_rate_cond_round_c")
nb_align_rate_cond_round_log = brm(num_gestural_align | rate(num_iconic_gestures) ~
1 + condition * log_round_c + num_lex_align +
(1+log_round_c|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope_gest_align_prop,
data = df_gest_align_posreg_prop,
sample_prior = T,
save_pars = save_pars(all = TRUE),
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/nb_align_rate_cond_round_log")
### loo compare
if (!file.exists("models/loo_nb_align_rate_cond_round.rds")){
nb_cond_round_loo = loo(nb_align_rate_cond_round)
saveRDS(nb_cond_round_loo, file = "models/loo_nb_align_rate_cond_round.rds")
}
if (!file.exists("models/loo_nb_align_rate_cond_round_c.rds")){
nb_cond_round_c_loo = loo(nb_align_rate_cond_round_c)
saveRDS(nb_cond_round_c_loo, file = "models/loo_nb_align_rate_cond_round_c.rds")
}
if (!file.exists("models/loo_nb_align_rate_cond_round_log.rds")){
nb_cond_round_log_loo = loo(nb_align_rate_cond_round_log)
saveRDS(nb_cond_round_log_loo, file = "models/loo_nb_align_rate_cond_round_log.rds")
}
nb_cond_round_loo = readRDS("models/loo_nb_align_rate_cond_round.rds")
nb_cond_round_c_loo = readRDS("models/loo_nb_align_rate_cond_round_c.rds")
nb_cond_round_log_loo = readRDS("models/loo_nb_align_rate_cond_round_log.rds")
loo_compare(nb_cond_round_loo, nb_cond_round_c_loo, nb_cond_round_log_loo)## elpd_diff se_diff
## nb_align_rate_cond_round 0.0 0.0
## nb_align_rate_cond_round_log -131.4 16.3
## nb_align_rate_cond_round_c -238.7 22.7The leave-one-out (LOO) Effect shows that the backward-difference coded round provides a far larger predictive power (elpd_diff) and a far smaller standard error (se_diff) compared to the centered round or centered log-transformed round. Thus, we will use the bd coded round as a predictor for further analyses.
8.4.1.2 Model Effect: ZI or not
zinb_align_rate_cond_round = brm(num_gestural_align ~
1 + condition * round + num_lex_align +
offset(log(num_iconic_gestures)) +
(1+round|pair) + (1|target),
family = zero_inflated_negbinomial(),
prior = priors_rslope_gest_align_prop_zinb,
data = df_gest_align_posreg_prop,
sample_prior = T,
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/zinb_align_rate_cond_round")
### loo compare
if (!file.exists("models/loo_zinb_align_rate_cond_round.rds")){
zinb_cond_round_c_loo = loo(zinb_align_rate_cond_round)
saveRDS(zinb_cond_round_c_loo, file = "models/loo_zinb_align_rate_cond_round.rds")
}
nb_cond_round_c_loo = readRDS("models/loo_nb_align_rate_cond_round.rds")
zinb_cond_round_c_loo = readRDS("models/loo_zinb_align_rate_cond_round.rds")
loo_compare(nb_cond_round_c_loo, zinb_cond_round_c_loo)## elpd_diff se_diff
## nb_align_rate_cond_round 0.0 0.0
## zinb_align_rate_cond_round -1.5 0.6The loo comparsion shows that non-inflation model has a higher predictive power than the zero-inflation model. As such, we will use NB regression models for further analyses.
8.4.1.3 Prior predictive check
nb_gest_align_prop_prior = brm(num_gestural_align | rate(num_iconic_gestures) ~
1 + condition + round + num_lex_align +
(1+round|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope_gest_align_prop,
data = df_gest_align_posreg_prop,
sample_prior = "only",
control = list(adapt_delta = 0.9,
max_treedepth = 20),
file = "models/nb_gest_align_prop_prior")
pp_check(nb_gest_align_prop_prior, ndraws = 100, type = "bars") +
coord_cartesian(xlim = c(0, 20),
ylim = c(0, 3000))
The prior predictive check shows that the model expects fewer amount of 2, 3, and 4. This suggests that the zero-inflation prior may be too large or the mean for the intercept prior is too low. We will check the prior-posterior update plot and posterior predictive check to see if the model generates data that are similar to the observed data. If not, we will consider modifying the priors.
8.4.1.4 Fit the model
nb_align_rate_cond_round = brm(num_gestural_align | rate(num_iconic_gestures) ~
1 + condition * round + num_lex_align +
(1+round|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope_gest_align_prop,
data = df_gest_align_posreg_prop,
sample_prior = T,
save_pars = save_pars(all = TRUE),
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/nb_align_rate_cond_round")
model = nb_align_rate_cond_round
summary(model)## Family: negbinomial
## Links: mu = log; shape = identity
## Formula: num_gestural_align | rate(num_iconic_gestures) ~ 1 + condition * round + num_lex_align + (1 + round | pair) + (1 | target)
## Data: df_gest_align_posreg_prop (Number of observations: 2199)
## Draws: 4 chains, each with iter = 20000; warmup = 2000; thin = 1;
## total post-warmup draws = 72000
##
## Multilevel Hyperparameters:
## ~pair (Number of levels: 45)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.67 0.09 0.51 0.87 1.00 20129
## sd(roundR12) 0.64 0.12 0.43 0.90 1.00 32271
## sd(roundR23) 0.13 0.08 0.01 0.30 1.00 20697
## sd(roundR34) 0.08 0.07 0.00 0.25 1.00 31671
## sd(roundR45) 0.13 0.09 0.01 0.34 1.00 25468
## sd(roundR56) 0.13 0.10 0.00 0.36 1.00 35078
## cor(Intercept,roundR12) 0.03 0.21 -0.38 0.43 1.00 42101
## cor(Intercept,roundR23) 0.18 0.32 -0.49 0.74 1.00 63141
## cor(roundR12,roundR23) -0.22 0.30 -0.73 0.44 1.00 56062
## cor(Intercept,roundR34) -0.00 0.33 -0.63 0.63 1.00 86753
## cor(roundR12,roundR34) -0.09 0.32 -0.67 0.56 1.00 75896
## cor(roundR23,roundR34) 0.01 0.33 -0.63 0.63 1.00 65940
## cor(Intercept,roundR45) 0.04 0.33 -0.60 0.65 1.00 79946
## cor(roundR12,roundR45) 0.09 0.31 -0.53 0.65 1.00 71978
## cor(roundR23,roundR45) 0.01 0.33 -0.62 0.64 1.00 57760
## cor(roundR34,roundR45) -0.04 0.34 -0.66 0.61 1.00 53393
## cor(Intercept,roundR56) -0.06 0.34 -0.68 0.60 1.00 85023
## cor(roundR12,roundR56) 0.05 0.32 -0.57 0.64 1.00 79246
## cor(roundR23,roundR56) 0.00 0.33 -0.63 0.63 1.00 65196
## cor(roundR34,roundR56) -0.01 0.33 -0.64 0.63 1.00 58872
## cor(roundR45,roundR56) -0.02 0.33 -0.65 0.62 1.00 56595
## Tail_ESS
## sd(Intercept) 38159
## sd(roundR12) 47486
## sd(roundR23) 26398
## sd(roundR34) 34620
## sd(roundR45) 31341
## sd(roundR56) 32642
## cor(Intercept,roundR12) 50684
## cor(Intercept,roundR23) 50618
## cor(roundR12,roundR23) 49539
## cor(Intercept,roundR34) 52609
## cor(roundR12,roundR34) 53428
## cor(roundR23,roundR34) 57498
## cor(Intercept,roundR45) 53014
## cor(roundR12,roundR45) 54864
## cor(roundR23,roundR45) 57810
## cor(roundR34,roundR45) 56422
## cor(Intercept,roundR56) 53281
## cor(roundR12,roundR56) 54104
## cor(roundR23,roundR56) 56138
## cor(roundR34,roundR56) 58815
## cor(roundR45,roundR56) 58976
##
## ~target (Number of levels: 16)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.06 0.04 0.00 0.14 1.00 22309 20923
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept -1.55 0.11 -1.77 -1.33 1.00 13026
## conditionAO_Asym 0.38 0.23 -0.08 0.84 1.00 16239
## conditionAsym_Sym 0.39 0.22 -0.06 0.83 1.00 12812
## roundR12 1.61 0.14 1.33 1.90 1.00 38765
## roundR23 0.24 0.08 0.09 0.39 1.00 54571
## roundR34 0.02 0.08 -0.14 0.19 1.00 58090
## roundR45 -0.12 0.10 -0.32 0.08 1.00 54198
## roundR56 -0.04 0.12 -0.27 0.19 1.00 64843
## num_lex_align -0.01 0.01 -0.04 0.02 1.00 100913
## conditionAO_Asym:roundR12 -0.07 0.28 -0.62 0.49 1.00 43120
## conditionAsym_Sym:roundR12 0.13 0.26 -0.39 0.64 1.00 36961
## conditionAO_Asym:roundR23 0.01 0.17 -0.32 0.34 1.00 57584
## conditionAsym_Sym:roundR23 -0.01 0.15 -0.30 0.29 1.00 55291
## conditionAO_Asym:roundR34 0.04 0.19 -0.33 0.40 1.00 55087
## conditionAsym_Sym:roundR34 0.15 0.16 -0.17 0.47 1.00 54834
## conditionAO_Asym:roundR45 0.24 0.22 -0.18 0.67 1.00 59488
## conditionAsym_Sym:roundR45 -0.14 0.19 -0.51 0.23 1.00 58188
## conditionAO_Asym:roundR56 0.19 0.26 -0.32 0.70 1.00 70505
## conditionAsym_Sym:roundR56 0.11 0.21 -0.30 0.51 1.00 72689
## Tail_ESS
## Intercept 26179
## conditionAO_Asym 30017
## conditionAsym_Sym 24995
## roundR12 48607
## roundR23 50032
## roundR34 50833
## roundR45 51064
## roundR56 54866
## num_lex_align 55847
## conditionAO_Asym:roundR12 52077
## conditionAsym_Sym:roundR12 47901
## conditionAO_Asym:roundR23 55257
## conditionAsym_Sym:roundR23 51500
## conditionAO_Asym:roundR34 55384
## conditionAsym_Sym:roundR34 54773
## conditionAO_Asym:roundR45 53403
## conditionAsym_Sym:roundR45 54752
## conditionAO_Asym:roundR56 55218
## conditionAsym_Sym:roundR56 57241
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 77.70 29.04 31.65 143.21 1.00 96781 51091
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).bayestestR::hdi(model)# bayestestR::hdi(model, ci = 0.89)The coefficients show that the proportion of gestural alignment was significantly higher in the SymAV condition than in the AO condition. Also, the rate of gestural alignment significantly increased from R1–R2 and R2–R3 and stabilized afterwards.
8.4.1.5 Hypothesis testing: Bayes factor
### varying priors for sensitivity analysis
prior_size = c("xs", "s", "l", "xl")
prior_sd = c(0.1, 0.3, 0.7, 1)
bfs_cond_ao_asym = c()
bfs_cond_asym_sym = c()
bfs_round12 = c()
bfs_round23 = c()
bfs_round34 = c()
bfs_round45 = c()
bfs_round56 = c()
bfs_lex_align = c()
for (i in 1:length(prior_sd)){
priors = c(
prior(normal(-1.39, 0.5), class = Intercept),
set_prior(paste0("normal(0,", prior_sd[i], ")"), class = "b"),
prior(normal(0, 0.5), class = sd),
prior(lkj(2), class = cor),
prior(normal(0, 50), class = shape))
fname = paste0("models/nb_align_rate_cond_round_", prior_size[i])
fit = brm(num_gestural_align | rate(num_iconic_gestures) ~
1 + condition * round + num_lex_align +
(1+round|pair) + (1|target),
family = negbinomial(),
prior = priors,
data = df_gest_align_posreg_prop,
sample_prior = T,
save_pars = save_pars(all = TRUE),
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = fname)
### BF for visibility conditions
# sym - asym
h = hypothesis(fit, "conditionAO_Asym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_cond_ao_asym = c(bfs_cond_ao_asym, bf)
# sym - ao
h = hypothesis(fit, "conditionAsym_Sym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_cond_asym_sym = c(bfs_cond_asym_sym, bf)
### BF for rounds
# R1 - R2
h = hypothesis(model, "roundR12 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round12 = c(bfs_round12, bf)
# R2 - R3
h = hypothesis(model, "roundR23 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round23 = c(bfs_round23, bf)
# R3 - R4
h = hypothesis(model, "roundR34 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round34 = c(bfs_round34, bf)
# R4 - R5
h = hypothesis(model, "roundR45 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round45 = c(bfs_round45, bf)
# R5 - R6
h = hypothesis(model, "roundR56 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round56 = c(bfs_round56, bf)
### BF for N. lex alignment
h = hypothesis(fit, "num_lex_align = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_lex_align = c(bfs_lex_align, bf)
}
### add BF for the main/medium model
prior_size[3:5] = c("m", prior_size[3:4])
prior_sd[3:5] = c(0.5, prior_sd[3:4])
### BF for visibility
# sym - asym
h = hypothesis(model, "conditionAO_Asym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_cond_ao_asym[3:5] = c(bf, bfs_cond_ao_asym[3:4])
# sym - ao
h = hypothesis(model, "conditionAsym_Sym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_cond_asym_sym[3:5] = c(bf, bfs_cond_asym_sym[3:4])
### BF for rounds
# R1 - R2
h = hypothesis(model, "roundR12 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round12[3:5] = c(bf, bfs_round12[3:4])
# R2 - R3
h = hypothesis(model, "roundR23 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round23[3:5] = c(bf, bfs_round23[3:4])
# R3 - R4
h = hypothesis(model, "roundR34 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round34[3:5] = c(bf, bfs_round34[3:4])
# R4 - R5
h = hypothesis(model, "roundR45 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round45[3:5] = c(bf, bfs_round45[3:4])
# R5 - R6
h = hypothesis(model, "roundR56 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round56[3:5] = c(bf, bfs_round56[3:4])
### BF for N. lex alignment
h = hypothesis(model, "num_lex_align = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_lex_align[3:5] = c(bf, bfs_lex_align[3:5])
### make a df for BFs
df_bf = data.frame(size = prior_size,
sd = prior_sd,
ao_asym = bfs_cond_ao_asym,
asym_sym = bfs_cond_asym_sym,
round12 = bfs_round12,
round23 = bfs_round23,
round34 = bfs_round34,
round45 = bfs_round45,
round56 = bfs_round56,
lex_align = bfs_lex_align) %>%
mutate(prior = paste0("N(0, ", sd, ")")) %>%
pivot_longer(cols = c("ao_asym", "asym_sym",
"round12", "round23", "round34", "round45", "round56",
"lex_align"),
names_to = "Effect",
values_to = "BF10") %>%
mutate(Predictor = ifelse(grepl("round", Effect), "Round",
ifelse(Effect == "lex_align", "N. lex align",
"Visibility")))
df_bf$Effect = recode(df_bf$Effect,
ao_asym = "AO--AsymAV",
asym_sym = "AsymAV--SymAV",
round12 = "R1--R2",
round23 = "R2--R3",
round34 = "R3--R4",
round45 = "R4--R5",
round56 = "R5--R6",
lex_align = "N. lex align")#### Plot BFs ####
ggplot(filter(df_bf, Effect!="R1--R2"), #exclude R1--R2 because BF is too huge
aes(x = factor(sd), y = BF10, group = Effect)) +
geom_hline(yintercept = 1, linetype="dashed") +
geom_point(aes(color=Effect)) +
geom_line(aes(color=Effect)) +
facet_wrap(vars(Predictor)) +
theme_clean(base_size = 15) +
theme(axis.text.x = element_text(colour = "black", size = 14),
axis.text.y = element_text(colour = "black", size = 14),
axis.title = element_text(size = 15, face = 'bold'),
axis.title.x = element_text(vjust = -2),
axis.title.y = element_text(vjust = 2),
legend.position = "top",
strip.text = element_text(size = 15, face = 'bold'),
plot.background = element_blank(),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines")) +
scale_y_log10("Bayes factor (BF10)",
breaks = c(0.001, 0.03, 0.01, 0.1, 0.33, 1, 3, 10, 30, 100),
labels = c(0.001, 0.03, 0.01, 0.1, 0.33, 1, 3, 10, 30, 100)) +
xlab("SD for the prior")
8.4.1.6 Probability of direction
p_direction(model)8.4.1.7 Pair-wise effect
emmeans(model, pairwise ~ condition)$contrasts## contrast estimate lower.HPD upper.HPD
## SymAV - AsymAV 0.388 -0.0568 0.829
## SymAV - AO 0.770 0.2787 1.265
## AsymAV - AO 0.382 -0.0787 0.832
##
## Results are averaged over the levels of: round
## Point estimate displayed: median
## Results are given on the log (not the response) scale.
## HPD interval probability: 0.95emmeans(model, pairwise ~ condition, level = 0.89)$contrasts## contrast estimate lower.HPD upper.HPD
## SymAV - AsymAV 0.388 0.028 0.744
## SymAV - AO 0.770 0.382 1.182
## AsymAV - AO 0.382 0.012 0.750
##
## Results are averaged over the levels of: round
## Point estimate displayed: median
## Results are given on the log (not the response) scale.
## HPD interval probability: 0.898.4.1.8 Prior-posterior update plot
post_sample = as_draws_df(model)
pp_update_plot(post_sample, model_type="nb")
8.4.1.9 Posterior predictive check
pp_check_overall = pp_check(model, ndraws = 100, type = "bars") +
coord_cartesian(xlim = c(0, 20),
ylim = c(0, 3000))
pp_check_sym = pp_check_each_condition(model, df_gest_align_posreg_prop, "SymAV")
pp_check_asym = pp_check_each_condition(model, df_gest_align_posreg_prop, "AsymAV")
pp_check_ao = pp_check_each_condition(model, df_gest_align_posreg_prop, "AO")
gridExtra::grid.arrange(pp_check_overall, pp_check_sym,
pp_check_asym, pp_check_ao,
ncol = 2)
Although the model prediction is not perfect, this model had a higher predictive power than the negative binomial model. As such, we will use this model.
8.4.1.10 Visualize the model estimates
plot(conditional_effects(model), ask=FALSE)
8.4.2 Model 6: [pair] condition_sum * round + num_lex_align
Here, we will run another model to compare the rate of gestural alignment between the SymAV and AO conditions.
8.4.2.1 Prior predictive check
nb_gest_align_prop_sum_prior = brm(num_gestural_align | rate(num_iconic_gestures) ~
1 + condition_sum + round + num_lex_align +
(1+round|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope_gest_align_prop,
data = df_gest_align_posreg_prop,
sample_prior = "only",
control = list(adapt_delta = 0.9,
max_treedepth = 20),
file = "models/nb_gest_align_prop_sum_prior")
pp_check(nb_gest_align_prop_sum_prior, ndraws = 100, type = "bars") +
coord_cartesian(xlim = c(0, 20),
ylim = c(0, 3000))
The prior predictive check shows that the model expects fewer amount of 2, 3, and 4. This suggests that the zero-inflation prior may be too large or the mean for the intercept prior is too low. We will check the prior-posterior update plot and posterior predictive check to see if the model generates data that are similar to the observed data. If not, we will consider modifying the priors.
8.4.2.2 Fit the model
nb_align_rate_cond_round_sum = brm(num_gestural_align | rate(num_iconic_gestures) ~
1 + condition_sum * round + num_lex_align +
(1+round|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope_gest_align_prop,
data = df_gest_align_posreg_prop,
sample_prior = T,
save_pars = save_pars(all = TRUE),
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/nb_align_rate_cond_round_sum")
model = nb_align_rate_cond_round_sum
summary(model)## Family: negbinomial
## Links: mu = log; shape = identity
## Formula: num_gestural_align | rate(num_iconic_gestures) ~ 1 + condition_sum * round + num_lex_align + (1 + round | pair) + (1 | target)
## Data: df_gest_align_posreg_prop (Number of observations: 2199)
## Draws: 4 chains, each with iter = 20000; warmup = 2000; thin = 1;
## total post-warmup draws = 72000
##
## Multilevel Hyperparameters:
## ~pair (Number of levels: 45)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.67 0.09 0.51 0.87 1.00 18575
## sd(roundR12) 0.64 0.12 0.43 0.90 1.00 31125
## sd(roundR23) 0.12 0.08 0.01 0.30 1.00 20919
## sd(roundR34) 0.08 0.07 0.00 0.24 1.00 32022
## sd(roundR45) 0.12 0.09 0.00 0.33 1.00 26026
## sd(roundR56) 0.13 0.10 0.00 0.37 1.00 34132
## cor(Intercept,roundR12) 0.04 0.21 -0.38 0.44 1.00 40280
## cor(Intercept,roundR23) 0.18 0.32 -0.49 0.73 1.00 65536
## cor(roundR12,roundR23) -0.22 0.31 -0.73 0.44 1.00 59711
## cor(Intercept,roundR34) -0.01 0.34 -0.64 0.63 1.00 93535
## cor(roundR12,roundR34) -0.09 0.32 -0.67 0.56 1.00 78483
## cor(roundR23,roundR34) 0.01 0.33 -0.62 0.64 1.00 67466
## cor(Intercept,roundR45) 0.04 0.33 -0.60 0.65 1.00 80878
## cor(roundR12,roundR45) 0.09 0.31 -0.54 0.66 1.00 78177
## cor(roundR23,roundR45) 0.01 0.33 -0.62 0.64 1.00 57583
## cor(roundR34,roundR45) -0.04 0.33 -0.66 0.61 1.00 54005
## cor(Intercept,roundR56) -0.06 0.34 -0.68 0.59 1.00 93756
## cor(roundR12,roundR56) 0.05 0.32 -0.57 0.64 1.00 85093
## cor(roundR23,roundR56) -0.00 0.33 -0.63 0.63 1.00 67602
## cor(roundR34,roundR56) -0.01 0.33 -0.63 0.62 1.00 60265
## cor(roundR45,roundR56) -0.02 0.33 -0.64 0.62 1.00 58058
## Tail_ESS
## sd(Intercept) 34700
## sd(roundR12) 48295
## sd(roundR23) 26610
## sd(roundR34) 31291
## sd(roundR45) 32316
## sd(roundR56) 33988
## cor(Intercept,roundR12) 50279
## cor(Intercept,roundR23) 51923
## cor(roundR12,roundR23) 48590
## cor(Intercept,roundR34) 54422
## cor(roundR12,roundR34) 54685
## cor(roundR23,roundR34) 57533
## cor(Intercept,roundR45) 51517
## cor(roundR12,roundR45) 55396
## cor(roundR23,roundR45) 57358
## cor(roundR34,roundR45) 57922
## cor(Intercept,roundR56) 54259
## cor(roundR12,roundR56) 55731
## cor(roundR23,roundR56) 56162
## cor(roundR34,roundR56) 59440
## cor(roundR45,roundR56) 59143
##
## ~target (Number of levels: 16)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.06 0.03 0.00 0.14 1.00 21296 20498
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept -1.54 0.11 -1.77 -1.33 1.00
## condition_sumAO_Sym 0.65 0.23 0.19 1.11 1.00
## condition_sumAsym_Sym 0.27 0.23 -0.18 0.72 1.00
## roundR12 1.62 0.14 1.33 1.90 1.00
## roundR23 0.24 0.08 0.09 0.39 1.00
## roundR34 0.03 0.08 -0.13 0.19 1.00
## roundR45 -0.12 0.10 -0.32 0.08 1.00
## roundR56 -0.03 0.12 -0.27 0.20 1.00
## num_lex_align -0.01 0.01 -0.04 0.02 1.00
## condition_sumAO_Sym:roundR12 0.03 0.28 -0.52 0.57 1.00
## condition_sumAsym_Sym:roundR12 0.13 0.26 -0.39 0.65 1.00
## condition_sumAO_Sym:roundR23 0.01 0.17 -0.32 0.34 1.00
## condition_sumAsym_Sym:roundR23 -0.01 0.15 -0.30 0.29 1.00
## condition_sumAO_Sym:roundR34 0.16 0.18 -0.20 0.52 1.00
## condition_sumAsym_Sym:roundR34 0.15 0.16 -0.17 0.47 1.00
## condition_sumAO_Sym:roundR45 0.13 0.22 -0.29 0.56 1.00
## condition_sumAsym_Sym:roundR45 -0.15 0.19 -0.52 0.22 1.00
## condition_sumAO_Sym:roundR56 0.25 0.26 -0.26 0.76 1.00
## condition_sumAsym_Sym:roundR56 0.08 0.21 -0.33 0.49 1.00
## Bulk_ESS Tail_ESS
## Intercept 13062 24757
## condition_sumAO_Sym 13128 25189
## condition_sumAsym_Sym 11201 23088
## roundR12 38899 48092
## roundR23 54168 50887
## roundR34 61182 52654
## roundR45 58074 53364
## roundR56 66159 55148
## num_lex_align 95708 56583
## condition_sumAO_Sym:roundR12 43962 53222
## condition_sumAsym_Sym:roundR12 36727 47437
## condition_sumAO_Sym:roundR23 58050 53368
## condition_sumAsym_Sym:roundR23 60204 55326
## condition_sumAO_Sym:roundR34 60163 56713
## condition_sumAsym_Sym:roundR34 60303 55921
## condition_sumAO_Sym:roundR45 60850 56172
## condition_sumAsym_Sym:roundR45 60868 55720
## condition_sumAO_Sym:roundR56 70417 57299
## condition_sumAsym_Sym:roundR56 74686 56413
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 77.56 29.03 31.62 143.29 1.00 107322 49549
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).bayestestR::hdi(model)# bayestestR::hdi(model, ci = 0.89)The coefficients show that the proportion of gestural alignment was significantly higher in the SymAV condition than in the AO condition. Also, the rate of gestural alignment significantly increased from R1–R2 and R2–R3 and stabilized afterwards.
8.4.2.3 Visualize posterior distributions
main_model = nb_align_rate_cond_round
model1 = zinb_align_cond_round
plot_posterior(main_model, model1, model)
8.4.2.4 Hypothesis testing: Bayes factor
### varying priors for sensitivity analysis
prior_size = c("xs", "s", "l", "xl")
prior_sd = c(0.1, 0.3, 0.7, 1)
bfs_cond_ao_sym = c()
for (i in 1:length(prior_sd)){
priors = c(
prior(normal(-1.39, 0.5), class = Intercept),
set_prior(paste0("normal(0,", prior_sd[i], ")"), class = "b"),
prior(normal(0, 0.5), class = sd),
prior(lkj(2), class = cor),
prior(normal(0, 50), class = shape))
fname = paste0("models/nb_align_rate_cond_round_sum_", prior_size[i])
fit = brm(num_gestural_align | rate(num_iconic_gestures) ~
1 + condition_sum * round + num_lex_align +
(1+round|pair) + (1|target),
family = negbinomial(),
prior = priors,
data = df_gest_align_posreg_prop,
sample_prior = T,
save_pars = save_pars(all = TRUE),
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = fname)
### BF for visibility conditions
# BF for sym - ao
h = hypothesis(fit, "condition_sumAO_Sym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_cond_ao_sym = c(bfs_cond_ao_sym, bf)
}
### add BF for the main/medium model
prior_size[3:5] = c("m", prior_size[3:4])
prior_sd[3:5] = c(0.5, prior_sd[3:4])
### BF for visibility conditions
# sym - ao
h = hypothesis(model, "condition_sumAO_Sym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_cond_ao_sym[3:5] = c(bf, bfs_cond_ao_sym[3:4])
### make a df for BFs
df_bf_temp = data.frame(size = prior_size,
sd = prior_sd,
Effect = "AO--SymAV",
BF10 = bfs_cond_ao_sym) %>%
mutate(prior = paste0("N(0, ", sd, ")"),
Predictor = "Visibility")
df_bf_new = df_bf %>%
filter(Effect != "N. lex align") %>%
rbind(df_bf_temp) %>%
rbind(df_bf_lex) %>%
mutate(Effect = factor(Effect,
levels = c("AO--SymAV", "AO--AsymAV", "AsymAV--SymAV",
"R1--R2", "R2--R3", "R3--R4", "R4--R5", "R5--R6",
"N. lex align")),
Predictor = factor(Predictor,
levels = c("Visibility", "Round", "N. lex align")))
df_bf_new %>% arrange(Effect, sd)#### Plot BFs ####
ggplot(filter(df_bf_new, Effect!="R1--R2"), #exclude R1--R2 because BF is too huge
aes(x = factor(sd), y = BF10, group = Effect)) +
geom_hline(yintercept = 1, linetype="dashed") +
geom_point(aes(color=Effect)) +
geom_line(aes(color=Effect)) +
facet_wrap(vars(Predictor), scales = "free_x") +
theme_clean(base_size = 15) +
theme(axis.text.x = element_text(colour = "black", size = 14),
axis.text.y = element_text(colour = "black", size = 14),
axis.title = element_text(size = 15, face = 'bold'),
axis.title.x = element_text(vjust = -2),
axis.title.y = element_text(vjust = 2),
legend.position = "top",
legend.title=element_text(size=14),
legend.text=element_text(size=13),
strip.text = element_text(size = 15, face = 'bold'),
plot.background = element_blank(),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines")) +
scale_y_log10("Bayes factor (BF10)",
breaks = c(0.001, 0.03, 0.01, 0.1, 0.33, 1, 3, 10, 30, 100),
labels = c(0.001, 0.03, 0.01, 0.1, 0.33, 1, 3, 10, 30, 100)) +
xlab("SD for the prior")
8.4.2.5 Probability of direction
p_direction(model)8.4.2.6 Prior-posterior update plot
post_sample = as_draws_df(model)
pp_update_plot(post_sample, model_type="nb")
8.4.2.7 Posterior predictive check
pp_check_overall = pp_check(model, ndraws = 100, type = "bars") +
coord_cartesian(xlim = c(0, 20),
ylim = c(0, 3000))
pp_check_sym = pp_check_each_condition(model, df_gest_align_posreg_prop, "SymAV")
pp_check_asym = pp_check_each_condition(model, df_gest_align_posreg_prop, "AsymAV")
pp_check_ao = pp_check_each_condition(model, df_gest_align_posreg_prop, "AO")
gridExtra::grid.arrange(pp_check_overall, pp_check_sym,
pp_check_asym, pp_check_ao,
ncol = 2)
Although the model prediction is not perfect, this model had a higher predictive power than the negative binomial model. As such, we will use this model.
8.4.2.8 Visualize the model estimates
plot(conditional_effects(model), ask=FALSE)
9 =====Correlation btw lexical and gestural alignment=====
cor = pcor.test(df_trial_info$num_lex_align, df_trial_info$num_gestural_align, df_trial_info$log_round_c)
cor10 =====Session info=====
sessionInfo()## R version 4.2.2 (2022-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 22631)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] rstan_2.32.6 StanHeaders_2.32.7 ppcor_1.1 MASS_7.3-58.1
## [5] doParallel_1.0.17 iterators_1.0.14 foreach_1.5.2 emmeans_1.10.7
## [9] tidybayes_3.0.7 bayestestR_0.15.2 brms_2.22.0 Rcpp_1.0.12
## [13] plotly_4.10.4 rsvg_2.6.2 DiagrammeRsvg_0.1 DiagrammeR_1.0.11
## [17] dagitty_0.3-4 ggh4x_0.3.1 ggthemes_5.1.0 hypr_0.2.8
## [21] plotrix_3.8-4 lubridate_1.9.3 forcats_1.0.0 stringr_1.5.1
## [25] dplyr_1.1.4 purrr_1.0.2 readr_2.1.5 tidyr_1.3.1
## [29] tibble_3.2.1 ggplot2_3.5.2 tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] colorspace_2.1-0 estimability_1.5.1 QuickJSR_1.1.3
## [4] rstudioapi_0.17.1 farver_2.1.1 svUnit_1.0.6
## [7] bit64_4.0.5 mvtnorm_1.2-4 bridgesampling_1.1-2
## [10] codetools_0.2-18 cachem_1.0.8 knitr_1.49
## [13] bayesplot_1.11.1 jsonlite_1.8.8 ggdist_3.3.2
## [16] compiler_4.2.2 httr_1.4.7 backports_1.4.1
## [19] Matrix_1.5-1 fastmap_1.1.1 lazyeval_0.2.2
## [22] cli_3.6.2 visNetwork_2.1.2 htmltools_0.5.8.1
## [25] tools_4.2.2 coda_0.19-4.1 gtable_0.3.6
## [28] glue_1.7.0 reshape2_1.4.4 posterior_1.6.1
## [31] V8_6.0.6 jquerylib_0.1.4 vctrs_0.6.5
## [34] nlme_3.1-160 crosstalk_1.2.1 insight_1.1.0
## [37] tensorA_0.36.2.1 xfun_0.53 ps_1.7.6
## [40] timechange_0.3.0 lifecycle_1.0.4 scales_1.3.0
## [43] vroom_1.6.5 ragg_1.3.0 hms_1.1.3
## [46] Brobdingnag_1.2-9 inline_0.3.21 RColorBrewer_1.1-3
## [49] yaml_2.3.8 curl_6.2.1 gridExtra_2.3
## [52] loo_2.8.0 sass_0.4.9 stringi_1.8.3
## [55] checkmate_2.3.1 pkgbuild_1.4.6 boot_1.3-28
## [58] cmdstanr_0.8.1 rlang_1.1.3 pkgconfig_2.0.3
## [61] systemfonts_1.0.6 matrixStats_1.3.0 distributional_0.5.0
## [64] pracma_2.4.4 evaluate_1.0.3 lattice_0.20-45
## [67] rstantools_2.4.0 htmlwidgets_1.6.4 labeling_0.4.3
## [70] processx_3.8.4 bit_4.0.5 tidyselect_1.2.1
## [73] plyr_1.8.9 magrittr_2.0.3 R6_2.6.1
## [76] generics_0.1.3 pillar_1.10.1 withr_3.0.2
## [79] datawizard_1.0.1 abind_1.4-8 crayon_1.5.3
## [82] arrayhelpers_1.1-0 tzdb_0.4.0 rmarkdown_2.29
## [85] grid_4.2.2 data.table_1.15.4 digest_0.6.35
## [88] xtable_1.8-4 textshaping_0.3.7 stats4_4.2.2
## [91] RcppParallel_5.1.7 munsell_0.5.1 viridisLite_0.4.2
## [94] bslib_0.9.0